Correlation of CMB with large-scale structure: I. ISW Tomography and Cosmological 

Implications 
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We cross-correlate large scale structure (LSS) observations from a number of surveys with cos- 
mic microwave background (CMB) anisotropies from the Wilkinson Microwave Anisotropy Probe 
(WMAP) to investigate the Integrated Sachs- Wolfe (ISW) effect as a function of redshift, covering 
2 ~ 0.1 — 2.5. Our main goal is to go beyond reporting detections towards developing a reliable like- 
lihood analysis that allows one to determine cosmological constraints from ISW observations. With 
this in mind we spend a considerable amount of effort in determining the redshift-dependent bias and 
redshift distribution (6(z) x dN/dz) of these samples by matching with spectroscopic observations 
where available, and analyzing auto-power spectra and cross-power spectra between the samples. 
Due to wide redshift distributions of some of the data sets we do not assume a constant bias model, 
in contrast to previous work on this subject. We only use the LSS data sets for which we can extract 
such information reliably and as a result the data sets we use are 2-Micron All Sky Survey (2MASS) 
samples, Sloan Digital Sky Survey (SDSS) photometric Luminous Red Galaxies, SDSS photometric 
quasars and NRAO VLA Sky Survey (NVSS) radio sources. We make a joint analysis of all samples 
constructing a full covariance matrix, which we subsequently use for cosmological parameter fitting. 
We report a 3.7(t detection of ISW combining all the datasets. We do not find significant evidence for 
an ISW signal at z > 1, in agreement with theoretical expectation in ACDM model. We combine the 
ISW likelihood function with weak lensing of CMB (hereafter Paper II [ij) and CMB power spectrum 
to constrain the equation of state of dark energy and the curvature of the Universe. While ISW does 
not significantly improve the constraints in the simplest 6-parameter fiat ACDM model, it improves 
constraints on 7-parameter models with curvature by a factor of 3.2 (relative to WMAP alone) to 



-0.004 



+0.014 



and with dark energy equation of state by 15% to w = — 1.01_q'4q [posterior me- 
dian with "1(t" (16th-84th percentile) range]. A software package for calculating the ISW likelihood 
function can be downloaded at http : //www. astro .princeton. edu/""shirley/ISW_WL .html, 

PACS numbers: 98.80.Es, 95.36. +x, 98.65.Dx. 



I. INTRODUCTION 



The Cosmic Microwave Background (CMB) has pro- 
vided us with a wealth of cosmological information. The 
large-scale anisotropies were first discovered by the Dif- 
ferential Microwave Radiometer (DMR) on Cosmic Back- 
ground Explorer (COBE) satellite and the smaller- 
scale CMB anisotropies were subsequently measured by 
various ground-based/balloon-borne experiments. More 
recently, the Wilkinson Microwave Anisotropy Probe 
(WMAP) satellite 0, Q produced a cosmic variance lim- 
ited map of CMB anisotropies down to I ~ 400. The 
structure of the angular power spectrum when combined 
with other cosmological probes (such as Sloan Digital Sky 
Survey, Hubble Key Project and 2dF Galaxy Red- 
shift Survey 0), allows extremely precise measurements 
of the cosmological parameters of the ACDM model. 
While most of the fluctuations seen by WMAP and other 
CMB experiments were generated at the last surface of 
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scattering, structures formed at low redshift also leave 
imprints on the CMB. These anisotropies, such as the 
thermal Sunyaev-Zeldovich (tSZ) Q and kinetic Sunyaev 
Zeldovich effects (kSZ) [§], the Integrated Sachs- Wolfe 
(ISW) effect and gravitational lensing, contribute 

only slightly to the CMB power spectrum on scales mea- 
sured by WMAP, but they can be detected by cross- 
correlating the CMB with suitable tracers of the large 
scale structure. 

This is the first of two papers that measure the Inte- 
grated Sachs- Wolfe effect and gravitational lensing (Pa- 
per II) in cross-correlation. In this paper, we focus 
on large scale galaxy-temperature correlations and their 
large scale cosmological source, the Integrated Sachs- 
Wolfe (ISW) effect. The ISW effect results from the red- 
or blue-shifting of the CMB photons as they propagate 
through gravitational potential wells. As the potential 
wells of the Universe (i.e., the spatial metric) evolve, the 
energy gained by photons falling into the potential well 
does not cancel out the energy loss as photons climb out 
of the well. This is important at late times when the 
Universe is not matter dominated and the gravitational 
potential is time dependent. It is only significant on large 
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scales, since on small scales the amount of time spent by 
the photon in each coherence region of the gravitational 
potential is small and any small scale fluctuations will be 
smoothed out as the photon go through numerous poten- 
tial wells along the line of sight. 

To measure the above effect, we cross-correlate the 
CMB temperature anisotropics with maps of galaxies 
from the Two Micron All Sky Survey (2MASS), lumi- 
nous red galaxies (LRGs) and quasars from the Sloan 
Digital Sky Survey, and radio sources from the NRAO 
VLA Sky Survey (NVSS). This incorporates most of the 
LSS tracers used by previous efforts (ill , [l^ . [isl . [li . [lil . 
M, E M, H H M, M, M, Hi, [Hi to detect the ISW 
effect. Our goal in this work extends this previous liter- 
ature by going beyond detecting the ISW effect to mea- 
suring its redshift evolution and using that to constrain 
different cosmological models (e.g. the ISW effect due to 
spatial curvature occurs at significantly higher redshifts 
than that due to a cosmological constant). We therefore 
require a large redshift range (z ~ to 2.5) but with suf- 
ficient redshift resolution to unambiguously discern any 
redshift evolution of the signal. In addition, to draw ro- 
bust cosmological conclusions from an observed redshift 
evolution, we must constrain both the redshift distribu- 
tion and evolution of the bias with redshift for each of 
the samples; the simple assumption of constant bias is 
in most cases no longer sufficient. These considerations 
drive our survey selections; we discuss these in more de- 
tail in Sec. IVIIII Our final product is a likelihood code 
that can be applied to any cosmological model. In addi- 
tion to providing complementary constraints on standard 
cosmological parameters, we expect it can be a strong 
discriminator of the modified gravity models, which have 
very distinctive ISW predictions [2g |. 

We review the theory behind the ISW effect in SecUll 
The CMB and LSS data sets used are described in 
Sec. IIII[ the results of cross-correlating the two are in 
Sec. IIVI Sec. |V] and [Vl] constrain the redshift distri- 
butions of the samples, and possible systematic contam- 
ination of the cross-correlations. Sec. IVIII presents the 
cosmological implications of these results, and Sec. IVIIII 
summarizes our conclusions. The companion paper (Pa- 
per II) uses the same data sets to detect the weak lensing 
of the CMB. All of the theoretical predictions are made 
with WMAP 3 year parameters (Ob/i2=0.0223, ^ch'^ = 
0.128, = 0, h = 0.732, erg = 0.761) except in Sec- 
tion [V] or otherwise stated. 



II. 



THEORY 



We briefly review the ISW effect and its cross- 
correlation with the galaxy density (see also Refs. (l3 . 
[27l . [Hi). The temperature anisotropy due to the ISW 
effect is expressed as an integral of the time derivative of 
the gravitational potential over conformal time 77, 



Arisw(^) = 2 



(1) 



where r]r and r]Q are the conformal time at recombina- 
tion and today, respectively, and we ignored the effect of 
Thomson scattering suppression, which is negligible for 
the redshift range of interest here. For scales sufficiently 
within the horizon, the gravitational potential (j> is re- 
lated to the mass fluctuation S = 5p/p in Fourier space 
by the Poisson equation: 



(k,z) = -^:Ja„(l + z) 
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(2) 



where flm is the ratio of the matter density to the critical 
density today, Hq is the Hubble constant today, c is the 
speed of light, z is the redshift, and k is the comoving 
wave number. On large scales where the mass fluctuation 
S -^1, the perturbations grow according to linear theory 
S{k,z) = 5lk,0)D{z)/D{Q). 

We are interested in cross-correlating the temperature 
anisotropies, St, with the observed projected galaxy over- 
density g. The intrinsic angular galaxy fluctuations are 
given by: 
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where b{z) is an assumed scale-independent bias factor 
relating the galaxy overdensity to the mass overdensity, 
i.e. 6g = bS, n(z) is the normalized selection function, 
and x(z) is the comoving distance to redshift z. We fo- 
cus on the cross-spectrum of the galaxies with the CMB 
temperature fluctuation: 
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where P{k) is the matter power spectrum today as a 
function of the wave number fc, and the functions [g]i 
and \T]i are 
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and 
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The Limber approximation, which is quite accurate when 
£ is not too small {£ > 10), can be obtained from Eq. (0]) 
by setting P{k) = P{k = {i + l/2)/x{z)) and using 
the asymptotic formula that {2/tt) J k'^dkji{kx)ji{kx') = 
{^lx )^{x - X') (when ^ > 1). We flnd that the substi- 
tution k = {£ + 1/2) /x{z) is a better approximation to 
the exact expressions than k = i/x{z). This gives 
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The above discussion ignores the effects of gravita- 
tional lensing, which alters the expected signal through 
two competing effects - changing the flux limit of the 
survey as well as the observed galaxy density. Both of 
these effects can be thought of as altering the redshift 
distribution of the tracers, and so we defer the discussion 
to Sec. El 



III. DATA 

We describe the CMB and galaxy data sets used in our 
analysis below; these are summarized in Table [B The 
data sets not used in this paper are discussed further 
in the Sec. IVIIIj where we provide detailed explanations 
for the choices made. All large scale structure data were 
pixelized in the HEALPix system with the resolution and 
sky coverage shown in Table [B 



A. CMB temperature from WMAP 

The WMAP mission 0, 3 measured the all-sky maps 
of the CMB at multipoles up to £ ~ several hundred. We 
use the second public data release of the WMAP data 
with the first three years of observations. The all-sky 
CMB maps are constructed in the following bands: K 
(23 GHz), Ka (33 GHz), Q (41 GHz), V (61 GHz) and W 
(94 GHz). These maps are pixelized in the HEALPix ^ 
resolution 9 format with 3 145 728 pixels, each 47.2 sq. ar- 
cmin in area. These maps are not beam-deconvolved and 
this, with the scan strategy of WMAP, results in nearly 
uncorrelated Gaussian uncertainties on the temperature 
in each pixel We limit our analysis to Ka through 
W band as the K-band is heavily contaminated by the 
Galactic emission. We trim all masks with the WMAP 
KpO mask and point source mask to remove regions con- 
taminated by Galactic emission and point sources, leav- 
ing 76.8% (2 414 613 resolution 9 HEALPix pixels) of the 
sky for the ISW analysis. We choose not to use either the 
WMAP "Internal Linear Combination" (ILC) map or the 
foreground cleaned map to avoid a number of practical 
difficulties as these maps lose frequency dependence of 
the original maps and have complicated pixel-pixel noise 
correlations. 



B. Two Micron All Sky Survey (2MASS) 

We use galaxies from the Two Micron All Sky Survey 
(2MASS) Extended Source Catalog (XSC) [13, [U, [32| 
as mass tracers of the low redshift Universe. The me- 
dian redshift of these objects is ~ 0.1. We use K20, the 
ifs-band isophotal magnitude measured inside a circu- 
lar isophote with surface brightness of 20 mag arcsec"^, 
as our default flux measure. We extinction correct the 
magnitudes from the catalog using the reddening maps 




K20 = K20.raiu ^ Ak , (8) 

where Ak = 0.367£'(B - T^) [li|. Note that we ignore 
changes to the isophotal radius due to extinction. We 
remove regions with Ak > 0-05 in the dataset as the 
galaxy density starts to drop drastically. We visually 
inspects how the galaxy density changes with Ak and 
decide to cut with Ak > 0.05 as there is a drastic drop. 
There are 1 586 854 galaxies in the 2MASS XSC after 
removing known artifacts and sources in close proximity 
to a large galaxy {cc_flag 7^'a' and 'z') and requiring 
use_src ~ 1 (which rejects duplicate observations of the 
same part of the sky). The 2MASS XSC can miss objects 
near bright stars or overlapping artifacts, and so we used 
the XSC coverage map [30|| and masked out pixels with 
< 98% coverage, thus ^ 8% of the sky. 

We divided the 2MASS sample into 4 flux bins: 12.0 < 
K20 < 12.5, 12.5 < K20 < 13.0, 13.0 < K20 < 13.5, 
13.5 < K20 < 14.0. Note that the redshift distribution 
of these 4 bins actually overlap significantly. Our sam- 
ple selection for 2MASS is similar to Afshordi et al. (l3 | 
except the pixelization. 



C. Data from Sloan Digital Sky Survey (SDSS) 

The Sloan Digital Sky Survey has taken ugriz CCD im- 
ages of 10"* deg^ of the high-latitude sky [33|- A dedicated 
2.5m telescope [ll,!!^! at Apache Point Observatory im- 
ages the sky in photometric conditions [s^l in five bands 
(uqriz) [3^[39| using a drift-scanning, mosaic CCD cam- 
era [33|- All the data processing are done by completely 
automated pipelines, including astrometry, source iden- 
tification, photometry (40l . [4ll | , calibration [43. l43|. spec- 
troscopic target selection |44l. [isl. li^ . and spectroscopic 
fiber placement [i^l- The SDSS is well underway, and has 
produced seven major releases (isl. I49I. IHol. Isil. Isa. IH^. Is^l . 

In addition to constructing LRG and quasar maps, we 
constructed three additional maps that we use to reject 
region sheavily affected by poor seeing or stellar contam- 
ination. These include (i) a map of the full width at half- 
maximum (FWHM) of the point-spread function (PSF) 
in r band; (ii) a map of stellar density (18.0 < r < 18.5 
stars, smoothed with a 2 degree FMHM Gaussian); and 
(iii) a similar map using only the red stars {g — r > 1.4). 

All SDSS magnitudes used here are extinction- 
corrected using the maps of Ref. js^]. We use SDSS 
model magnitudes for the LRGs, and PSF magnitudes 
for the quasars and stars. 



1. Luminous Red Galaxies 

We use the photometric Luminous Red Galaxies 
(LRGs) from Sloan Digital Sky Survey (SDSS) con- 
structed as described in [55|. The LRGs have been very 
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TABLE I: The large-scale structure data sets used. The effective bias bcs and bias-weighted redshift {z)b are given here for 
the purpose of qualitatively illustrating which redshift ranges are probed by each sample. They are computed for the fiducial 
WMAP cosmology as b^s = / f{z) dz and {z)t = f zfjz) dz/h^ff, respectively; the redshift distributions f{z) will be computed 
in Sec.[Vl The data are pixelized using HEALPix |29l | at the resolutions listed in the table. 



Sample (its notation in paper) 


Area 


Density 


Number of 


HEALPix 


Number of 


bcff 






deg^ 


deg"^ 


galaxies 


resolution 


HEALPix Pixels 






2MASS, 12.0 <Ks< 12.5 (2MASS0) 


27191 


1.84 


50 096 


9 


2 073 457 


1.63 


0.06 


2MASS, 12.5 < 7^, < 13.0 (2MASS1) 


27191 


3.79 


103 060 


9 


2 073 457 


1.52 


0.07 


2MASS, 13.0 < Ks < 13.5 (2MASS2) 


27191 


7.85 


213 516 


9 


2 073 457 


1.54 


0.10 


2MASS, 13.5 <Ks < 14.0 (2MASS3) 


27191 


16.0 


435 570 


9 


2 073 457 


1.65 


0.12 


SDSS, LRG, low-z (LRGO) 


6 641 


35.1 


232 888 


10 


2 025 731 


1.97 


0.31 


SDSS, LRG, high-2 (LRGl) 


6 641 


93.8 


622 646 


10 


2 025 731 


1.98 


0.53 


SDSS, QSO, low-z (QSOO) 


6 039 


20.8 


125 407 


10 


1 842 044 


2.36 


1.29 


SDSS, QSO, high-2 (QSOl) 


6 039 


18.3 


110 528 


10 


1 842 044 


2.75 


1.67 


NVSS point sources (NVSS) 


27 361 


40.3 


1 104 983 


8 


521 594 


1.98 


1.43 



useful as a cosmological probe since they are typically 
the most luminous galaxies in the Universe, thus prob- 
ing a larger volume than most other tracers. On top of 
this, they also have very regular spectral energy distri- 
butions and a prominent 4000A break, making photo-z 
acquisition much easier than the other galaxies. We will 
not be repeat our selection criteria here as it is thor- 
oughly described in [HBl- We only accept sky regions 
with E{B -V) < 0.08 (almost identical to A < 0.2 as 
in [111) and an r band FWHM< 2.0 arcsec. 

Furthermore, there are a few regions in SDSS that 
have >60% more red stars than typical for their galactic 
latitude; we suspect photometric problems and rejected 
these regions. The red star cut removed 427 deg^ in as- 
sorted parts of the sky. 

We sHce our LRG sample into two redshift bins for the 
ISW analysis: 0.2 < Zphoto < 0.4 and 0.4 < Zphoto < 0.6. 



2. Photometric quasars 

We select quasars photometrically from the Sloan Dig- 
ital Sky Survey by first generatin g a candidate quasar 
catalog consisting of UVX objects [5a|. These are point 
sources with excess UV flux (i.e. u — g < 1.0) observed g 
magnitudes fainter than 14.5 (to avoid saturation prob- 
lems), extinction corrected g magnitudes brighter than 
21.0, and u-band error less than 0.5 mag (> 2a detection 
in u). We call this the ALL-UVX catalog. We also have 
the public catalog of photometric quasars from Data Re- 
lease 3 (DR3) generated by Ref. [53|, which we will call 
DR3-QS0 objects. We also construct a UVX object list 
from only DR3 data, denoted DR3-UVX. This catalog 
is used to extend the selection and photometric redshifts 
from the DR3 region to the ALL region. Ideally the cat- 
alog would have been based on running the algorithm 
of Ref. [131 on the ALL region but this option was not 
available at the time we constructed the quasar catalog. 

We first match the DR3-UVX objects to the DR3-QS0 
objects and then assign the photometric redshifts from 
the DR3-QS0 objects to the matched DR3-UVX object. 



For objects that are in DR3-UVX catalog, but not in the 
DR3-QS0 catalog, we mark them as rejects. We now 
have a DR3-UVX catalog with every object either as- 
signed a redshift or marked as a reject. The reject rate 
for DR3-UVX (ALL-UVX) is 89% (93%). Then, we lay 
down the DR3-UVX catalog in color^ {u—g,g—r,r—i,i~z) 
space, and then for each ALL-UVX object, we find its 
nearest neighbor in this color''^ space, then assigning it the 
same "redshift" as its matched DR3-UVX neighbor. If the 
DR3-UVX object has a redshift (not a reject), then the 
ALL-UVX object is classified as a quasar with the same 
redshift (photo-z only), otherwise it is rejected. This 
procedure generates a photometric catalog of quasars in 
the full survey area, based on the matching against DR3 
quasars in color* space. However, this catalog only has 
the photometric redshifts, but not the actual redshift dis- 
tribution. The actual redshift distribution will be dis- 
cussed in Sec.lVIl The average color offsets of the quasar 
candidate to its match for u — g, g — r, r~i and i ~ z 
are 0.0018, 0.0056, 0.0075 and 0.0045, while the typical 
errors on the colors of the candidates are 0.11 {u ~ g), 
0.13 {g - r), 0.14 (r - i) and 0.17 (j - z). As the color 
differences between the match and the candidate are well 
within the error of the colors, we conclude that the quasar 
candidates are matched with high accuracy. 

We then cut the catalog according to E{B — V)< 0.05 
and FWHM< 2.0 arcsec. These cuts are determined 
when we look at the variation of the quasar number 
overdensity over a range of extinction and seeing. Also, 
since quasars are more sensitive than LRGs to extinction 
(as a result of the importance of the u filter in select- 
ing quasars), we cut the catalog at a lower E{B — V). 
We also imposed a cut rejecting regions with more than 
twice average stellar density, i.e. we require rigtar < 564 
stars/deg^. 

We further divide the sample into two redshift (photo- 
z) bins: 0.65 < Zphoto < 1-45 (low-z) and 1.45 < Zphoto < 
2.0 (high-z). This division of sample is due to the fact 
that there are strong emission lines (e.g. Mg ii) that red- 
shift from one filter into the next around the redshifts 
of 0.65, 1.45 and 2.0, causing these two redshift bins to 



5 





FIG. 1: The overdensity maps of various tracer samples in Galactic coordinates. The scale runs from g = —1 (black, no 
galaxies) to g = —0.25 (blue), g = (green), g = +0.25 (red), and g = +1 (white, > 2x mean density). 



be relatively free of cross-contamination. However, as 
we will see, they do contain significant contamination 
from redshifts below 0.65 and above 2.0. We therefore 
constrain their redshift distribution by cross-correlating 
these with auxiliary data sets; we discuss this further in 
Sec. El 

The construction of the full sample using the DR3 cat- 
alog as described above introduces one potentially worry- 
ing systematic, namely the possibility that regions of the 
sky observed after DR3 would have a different density 
of sources than DR3 regions as a result of the nearest- 
neighbor method misbehaving in low-density regions of 
color^ space. This would provide a spurious feature in 
the quasar maps that resembles the DR3 coverage map. 
In order to check for this problem, we look for corre- 
lations between observing dates (if the ALL sample is 
misbehaving, it will be different from DR3 sample) with 
galaxy overdensity, and we do not find any significant 
correlations (Fig. [2]) . We also look at the correlation be- 
tween quasar overdensity and the stellar number density 
to see if there is significant stellar contamination, we do 
not find any either (Fig. [2|) . 



D. NRAO VLA Sky Survey (NVSS) 

The NRAO VLA Sky Survey (NVSS) is a lA GHz 
continuum survey covering the entire sky north of —40° 
declination using the compact D and DnC configurations 
of the Very Large Array (VLA) (H^]. The images all have 
45 arcsec FWHM resolution and nearly uniform sensitiv- 



ity and yield a catalog of almost 2 x 10^ discrete sources 
stronger than ^ 2 mJy. 

This survey has several potentially major artifacts: 
Galactic synchrotron emission, spurious power from 
bright sources and a declination-dependent striping prob- 
lem. All of these have to be treated properly before one 
can claim that the power coming from the cross-/auto- 
correlation is not due to some spurious issues. The Galac- 
tic synchrotron emission can in principle be an issue be- 
cause it contributes significantly to the noise temperature 
of the VLA, and for realistic number counts, increased 
noise temperature could change the number of sources 
with measured flux above some threshold. (As an inter- 
ferometer the VLA is not directly sensitive to the dif- 
fuse synchrotron foreground.) This issue is treated by 
incorporating a template - the Haslam map [H^ - in the 
cross-correlation analysis and projecting out the power 
that are correlated to this template. Even though the 
Haslam map is at 408 MHz, the frequency dependence 
of the galactic synchrotron emission is fairly flat, allow- 
ing us to use it as a template of the Galactic synchroton 
radiation. The bright sources are problematic since the 
VLA has a flnite dynamic range (~ 1000 in snapshot 
mode with limited itw-plane coverage) and thus the iden- 
tification of faint sources in fields with a bright source is 
unreliable. This issue is mitigated by masking out all the 
bright sources. Striping is a known systematic effect in 
NVSS [iOl: the galaxy density has a systematic depen- 
dence on declination, which can mimic long-wavelength 
modes in the galaxy field. To deal with the above po- 
tential problems, we first impose a flux limit of 2.5 mJy 
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FIG. 2: LRG and QSO overdensity vs various quantities such as reddening, PSF FWHM (r band), observing time (MJD), red 
star density, and star density. In each panel the circles show the low-redshift sample and the squares show the high-redshift 
sample. The Modified Julian Date (MJD) of the DR3 ending date is 52821. Note that there are very few accepted pixels at 
the extremes of reddening, PSF FWHM, and stellar density, resulting in the large fluctuations seen in the figure. 



(where NVSS is 50% complete), mask out a 0.6 degree 
radius around all the bright sources (> 2.5 Jy). Then 
to reduce striping, we also include templates to project 
out the synchrotron and declination-striping modes. The 
implementation of this projection of spurious power will 
be further discussed in Sec. IIVI 



IV. CROSS-CORRELATION POWER 
SPECTRUM ANALYSIS 

A. Methodology 

We start by organizing the temperature fluctuations 
and the galaxy overdensities into a single data vector. 



X = (XB^TjXg 



(9) 



where xb^t is a vector with the measured CMB temper- 
ature (with the monopole and dipole subtracted) in band 
B at every HEALPix pixel; analogously, Xg is the tracer 
number overdensity. The vector x has a total length 
A^pix,CMB+-/Vpix,LSS where iVpix.cMB and iVpi^.LSS are the 
number of accepted pixels for the CMB and LSS maps 
respectively. We suppress the band subscript for simplic- 
ity, with the implicit understanding that we always refer 
to the cross correlation of a single WMAP band with the 
tracer overdensity. The covariance matrix of x is, 



C = C 



diag 



Cf^t 





where Cdiag is given by. 



C99 + N99 



(10) 



(11) 



where N^^ is the noise matrix. 
C99 and C9'^ are deflned by 



The submatrices C 



TT 



i„{h1)YUn]) 



(12) 



where is the position (on the sky) of the i*^ point of the 
vector Xq. The temperature-temperature, galaxy-galaxy 
and galaxy-temperature angular power spectra are de- 
noted by Cf^^Cf^ and Cf'^ respectively. 

The galaxy power spectrum is flrst estimated using a 
pseudo-C/ estimator [6l||, and flt by the non- linear power 
spectrum of js^l, multiplied by a constant linear bias. We 
project out the monopole and dipole of both these power 
spectra by setting the power in the I = 0,1 modes to a 
value (10~^) much greater than the true power spectrum. 

We parametrize Cf as a sum of bandpowers, Pi,i, with 
amplitudes q to be estimated. 



1,1 



We consider "flat" bandpowers given by 

hi 



B{1) 




otherwise, 



(13) 



(14) 



where B(l) is the product of the beam transfer func- 
tion [g^I, and the HEALPix pixel transfer functions at 
WMAP and LSS resolution. This parametrizes the power 
spectrum as a sum of step functions and is useful when 
the shape of the power spectrum is unknown. 

We estimate the Ci by forming quadratic combinations 
of the data Hlfil, 



tag ^diag 



-1-1 



(15) 
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These are related to the estimated Ci by the response 
matrix F, 



where 



Fi-j = -tr 
y 2 



■ 1 dc , dc 

diag g^. diag g 



(16) 



(17) 



If Cf^ <C ^Jcf^Cf^ , then the are good approxima- 
tions to the maximum likelihood estimates of the The 
covariance matrix of the q is the inverse of the response 
matrix, if the fiducial power spectra and noise used to 
compute C^-^g correctly describe the data (in this case 
F is the Fisher matrix, hence the notation). The ma- 
trix Cdiag determines the weighting and is often called 
a "prior" in quadratic estimation theory. Note that this 
usage has nothing to do with Bayesian priors - in partic- 
ular, Eq. ifTB]) is unbiased regardless of the choice of prior 
(though for bad choices the estimator is not minimum 
variance). Implementing the above algorithm is compli- 
cated by the sizes of the datasets; the implementation we 
use is in [H, [IB, [H^l, and we refer to the reader to the 
discussion there. 

[In addition to the cross-power spectra in Eq. (fT4l) . 
in quadratic estimator theory one usually tries to esti- 
mate the CMB and galaxy auto-power spectra as well. 
Because our prior is diagonal, however, these decouple, 
i.e. the entries in that couple the auto-powers and 
cross-powers are zero. For this reason we can leave the 
auto-powers out of the quadratic estimator.] 

As mentioned earher, the NVSS dataset has issues that 
require additional processing. Assume a systematic E 
that we characterize as follows: 



B. Priors 

To generate the priors Cdiag for the cross-correlation 
power spectrum analysis, we need the approximate au- 
topower spectrum of the galaxies. The auto-correlation 
is done using the same methodology as described in 
Sec. IIV Al The resulting autopower spectra must be 
smoothed, before being used as priors. This avoids sta- 
tistical fluctuations in Ci over- or under-weighting the 
corresponding monopoles in the cross-correlation, which 
could result in underestimation of C^^ signal since we 
would artifically down-weight multipoles that had acci- 
dentally high power in galaxies and place more weight on 
multipoles that had little power. We did the smoothing 
in two different ways. For the cases where the redshift 
distribution was available early enough in the analysis 
(2MASS or LRG), we flt the auto-power spectrum to the 
non-linear matter power spectrum [g^l to get the Hnear 
bias. In other cases (quasars, NVSS) we did not have 
the redshift distribution at the time the priors were cre- 
ated; we created the priors by using a smoothed, splined 
auto-power spectrum of the sample as the prior. 

In the cases where we did a flt using the nonlinear mat- 
ter power spectrum, the fit biases are 1.15, 1.18, 1.20, 
and 1.22 (2MASS, brightest to faintest); 1.92 (LRG low- 
z)\ and 1.86 (LRG high-z). After generating the priors, 
we made several modifications to the analysis, includ- 
ing the inclusion of redshift-dependent bias in 2MASS. 
Thus while the priors were not updated since they give a 
good fit to the observed autopower spectrum, it should 
be noted that these bias values are not used in the cos- 
mological analysis (i.e. for ISW prediction purposes). 

To generate priors for the CMB, we generate the pri- 
ors using the theoretical C^s from WMAP and take into 
the account of the effect of pixelization and beams by 
convolving with the pixel and beam window functions. 



If estimate Ci, even if C is the true covariance, we will 
still have a biased answer. However, the substitution 

C ^ C*™ + (^EE^ (19) 

yields an unbiased estimate of q when ^ — > cx). One 
can add as many systematic templates E (i.e. modes 
to project out of the map) as desired. To immunize the 
NVSS correlations from possible systematics, we break 
the NVSS map into 74 decHnation rings, and for each 
ring include a template map E consisting of either -1-1 
(for pixels within the declination ring) or (for all other 
pixels). This removes the declination-dependent stripes. 
We also put in the 408 MHz Haslam map [5§| (techni- 
cally THasiam " 20 K) as a template for the Galactic syn- 
chrotron radiation. We experimented with the values of C, 
and found that the cross-spectra are converged with the 
choice C = 1 for the declination rings and C = 10~^ 
for the synchrotron map. 



C. Results of cross-correlation 

Figs. [31 m and [5] plot the cross-correlation between 
WMAP and the 2MASS, SDSS and NVSS samples re- 
spectively; the four different symbols in each of these 
plots correspond to the four WMAP bands we use. The 
observed achromatic nature of the signal is consistent 
with it being ISW, and is an important check for fre- 
quency dependent systematics. The two quasar samples 
are at the highest redshifts we can probe, so if there is 
an ISW cross-correlation at z ~ 1-2, it would mean that 
there is significant gravitational potential change at these 
redshifts. This is not expected in simplest ACDM cos- 
mology, but could be present either in models where dark 
energy equation of state is rapidly changing with redshift 
or in models where curvature plays a role. The observed 
lack of a signal for these redshifts therefore strongly con- 
strains such models. Note however that the NVSS cross- 
correlation cannot be automatically interpreted as a de- 
tection of high redshift ISW, as (see below) it covers a 
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FIG. 3: Galaxy density correlations with WMAP temperatures (4 bands: Ka (crosses), Q (triangles), V (squares), W (empty 
triangles), error bars are from the correlations with V-band. This contains 2MASS galaxy density correlations with WMAP, 
starting from (from left to right, top to bottom) the brightest sample, to the bottom the dimmest sample. We shift the points 
on X-axis for clarity. The dotted line shows the predicted signal for the sample with WMAP 3-year parameters and bdN/dz 
estimated in Sec. |Vl 



wide redshift range. 



V. REDSHIFT DISTRIBUTIONS 

The basic problem is to determine for each galaxy sam- 
ple i and each cosmological model the function fi{z) that 
relates the matter density S(r) to the two-dimensional 



galaxy overdensity g^: 

/>oo 

g,(n)= / Mz)S[n,x{z)]dz. (20) 
Jo 

Eq. l(20|) is understood to be valid on scales where the 
galaxies trace the matter distribution. In the absence 
of magnification bias, the function fi{z) is simply the 
product of the bias and the redshift distribution: fi{z) = 
bi{z)Ili(z) , where 11^(2;) is the probability distribution for 
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FIG. 4: Same as Fig. [3] except for the SDSS density maps from (from left to right, top to bottom): low-z LRG, high-z LRG, 
low-z QSO, high-z QSO. 



the galaxy redshift. In the presence of magnification bias, 
which is important for the SDSS quasars and possibly the 
NVSS radio sources, fi{z) takes on the more compHcated 
form 



Mz) = h{z)U,{z) + / W{z,z')[a{z') - l]U,{z')dz', 

J Z 

(21) 

where a{z') is the slope of the number counts of the 
galaxy density as a function of fiux: iV(> F) cx 



Here W{z, z') is the lensing window function: 

Wiz,z') - ^n„,H^^^^smlx{z) 

X [cotK x(z) - cotx x{z')] , (22) 

where xi^) = Jq dz" / H{z") is the radial comoving dis- 
tance, sinif X is the sine like function (equal to x in a flat 
Universe), and coti<- x = ^(Insinif x)/dx is the cotangent 
like function (equal to 1/x in a flat Universe) . 

It is in fact the function fi{z) that is required if one 
is to predict the ISW effect in a given cosmology. It is 
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FIG. 5: Same as Fig.[3]except for the NVSS cross-correlation. 



this same function that is required to predict the linear- 
regime angular power spectrum of the galaxies. This sec- 
tion describes the method by which fi{z) is obtained for 
each of the samples. The methods are quite different due 
to the different types of information available for each 
sample. In particular there are very few spectroscopic 
redshifts available for NVSS. Note however that all meth- 
ods include galaxy clustering data, as this is needed to 
determine the bias even if the redshift probability distri- 
bution Tli{z) is known perfectly. 

All of the numbers and plots in this section only that 
depend on cosmology are computed using the original 
WMAP third-year fiat 6-parameter ACDM cosmology 
{nbh^ = 0.0222, n„,h'^ = 0.1275, h = 0.727, erg = 0.743, 
and Us = 0.948), i.e. from the first release of Ref. [63 |. 
However in the Markov chain, the functon fi{z) is re- 
computed for each cosmological model and used to pre- 
dict the ISW signal. 



A. 2MASS 

The 2MASS samples go down to a limiting magnitude 
of K20 = 14. At this relatively bright magnitude, al- 
most all objects (97.9%, after correcting for the fiber col- 
lisions) have SDSS spectra, provided of course that they 
lie within the spectroscopic mask. In practice there are 
two subtleties that can occur. One is that the bias 52mass 
cannot be obtained to high accuracy from linear theory 
because even the moderate multipoles {I ~ 20) are non- 
linear, especially for the nearest 2MASS sHce, and the 
lowest multipoles suffer from cosmic variance. The other 
is that the bias varies with redshift: even though the 
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FIG. 6: The 2MASS redshift distribution, binned in units of 
Az — 0.01. The top panel shows the raw measured distribu- 
tion, n(z), and the bottom panel is corrected for relative bias 



2MASS galaxies cover a narrow range in redshift during 
which the Universe expands by only ^ 30%, the use of 
apparent magnitude to define the samples means that 
the typical luminosity of a galaxy varies by several mag- 
nitudes across the redshift range of interest, more biased, 
this effect shifts the peak of the effective redshift distribu- 
tion f{z) to higher redshifts than the actual distribution 

n(z). 

We match the 2MASS galaxies with the SDSS MAIN 
galaxy sample by first defining the 2MASS sample as 
discussed in IIIIBl then we select 2MASS galaxies only 
within mask that is more than 90% complete. We then 
try to match all the 2MASS galaxies with the SDSS 
MAIN galaxies that are within 3" and found that al- 
most all of the objects from 2MASS sample have SDSS 
spectra. We thus use the spectroscopic redshifts of the 
matched SDSS galaxies to identify the redshifts of the 
2MASS galaxies. The redshift distribution is binned with 
Sz — 0.01. The redshift distribution for each of the four 
slices is shown in Fig. [6l 

The problem of nonlinear evolution is generally very 
complicated, however for ISW work we only need a solu- 
tion accurate to a few tens of percent. Therefore we have 
used the Q-model Q, which relates the galaxy power 
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spectrum to the linear power spectrum via 



Galaxy power spectra 



,1 



1 + Ak 



(23) 



where b is the linear bias appearing in Eq. l(2T|) . Cole et al. 
Q found in simulations that this function fits the galaxy 
power spectrum in simulations for A = 1.7/i~^ Mpc, 
while the required value of Q varies depending on the 
sample. Our method is to compute the theoretical angu- 
lar galaxy power spectrum C^^(th) via the Limber inte- 
gral, and fit this to the measured C^^ treating b and Q 
as free parameters. This procedure can be done either 
assuming b is constant with redshift, or (better) taking 
into account the redshift-dependent bias, 



^'gal(fc,2) 



(^) 



1 + Qk^ 
1 + Ak 



Pii„(fc,z), (24) 



where b^ei{z) is known and 6* is a free parameter. While 
there is very little evolution in the 2MASS redshift range, 
the nearby and distant galaxies can have very different 
biases because they correspond to different luminosity 
ranges. The results for each are shown in Table Ull b^diz) 
is based on taking the r-band luminosities of the galaxies 
and using 6rci(i) from Tegmark et al. @. Note that 
the prominent peak of redshift distribution at z ~ 0.08 
is a supercluster known as the Sloan Great Wall. (In 
principle Q can depend on redshift as well, so one should 
be careful about interpreting the fit value and indeed one 
can see from Table [TTl that Q fit in this way is not stable. 
However the < Icr changes in (6) seen in the table when 
we restrict to much lower ^max suggest that this is not a 
large effect on the bias.) 

The Q-model fits for the 2MASS sample (and the 
LRGs) are shown in Fig. [71 



B. SDSS LRGs 

Next we consider the photometric LRG sample from 
SDSS. The sample is faint enough that spectroscopic 
redshifts are unavailable for most of the objects. For- 
tunately, precise photometric redshifts are available for 
LRGs since they have very uniform spectra whose main 
broadband feature is a break at 400 nm. This break 
passes through the SDSS g and r filters in the interest- 
ing redshift range, so the g ~ r and r — i colors of an 
LRG correlate very strongly with its redshift [53|. The 
error distribution of the photometric redshifts has been 
calibrated using spectro-zs from the 2SLAQ survey [g^; 
this procedure, and an inversion method used to deter- 
mine the actual redshift distribution given the photo-z 
distribution, are described in Padmanabhan et al. [55| . 
These methods were applied to determine the redshift 
probability distribution lli{z) for the LRGs used in this 
sample. The redshift distributions so obtained are shown 
in Fig. E 
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FIG. 7: The galaxy power spectra for the four 2MASS and 
two SDSS LRG samples, and the Q-model fits. The solid 
lines show the range of multipoles used in the fit, the dashed 
lines are extrapolations. Note that at very small scales the 
Q-model is not a good description of the power spectrum. 



The bias is determined by the same Q-model fitting 
procedure as we used for 2MASS. The maximum val- 
ues of £ considered are 240 for the low-z slice and 400 
for the high-z slice, which correspond to roughly k sa 
0.3/iMpc~^ at the typical redshifts of these samples. For 
the fiduciual cosmology, the low-z LRG slice gives a bias 
of 6 = 1.97 ± 0.05 and Q = 21.7 ± 2.6; the high-z slice 
gives b = 1.98 ± 0.03 and Q = 17.1 ± 1.5. In order to re- 
duce the possible impact of the nonhnear regime on our 
results, we also did fits where the maximum value of £ 
was reduced by a factor of 2 or 4. The results are shown 
in Table Hill and the bias estimates are seen to be consis- 
tent with each other. In what follows we have used the 
original (Zmax = 240,400) fits for the LRG bias, noting 
that the remaining uncertainty in b is small compared 
to the uncertainty (change in number of sigma detection 
is: 0.0043 (0.0388) for low-z LRG (high-z LRG)) result- 
ing from statistical error in the ISW signal. However we 
note that it is not clear how well the Q-model works for 
LRGs at small scales, and we recommend more detailed 
analysis before taking the very small statistical error in 
b at face value. The Q-model fits are shown in Fig. \7\ 

For the LRGs - unlike the 2MASS galaxies - each of 
the two photo-z slices covers a narrow redshift range and 
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TABLE II: The bias of the 2MASS galaxies as determined using the Q-model parametrization. The second column in each 
line shows the maximum value of I used in the main fits (varying or constant b). The first fit ("varying 6") uses Eq. l[24|l and 
should be viewed as the main result. For this fit we show the mean bias, i.e. {&) = / h{z)\i{z) dz, as this is easier to compare 
with other results than 6*. The second fit ("constant 6") has the bias fixed to a constant value. The third fit (/max = 24) has 
a bias varying according to Eq. I|24p but the fit is restricted to the region / < 25 in order to reduce the effect of nonlinearities. 
Note that the biases obtained from the varying-fc fits are consistent with each other, while the constant-bias fit finds a lower 
value of h by up to ~ 6% depending on the sample. 
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LRG redshift distributions 




FIG. 8: The redshift distributions of the two LRG samples. 
The dashed lines show the probability distribution for the 
photo-2s, where as the solid ("deconvolved") lines show the 
smoothed true redshift distribution based on the reconstruc- 
tion method of Padmanabhan et al. ISSll. 



TABLE III: The LRG bias and Q-parameter determined using 
several ranges of I. The "original" value /orig is 240 for the low- 
z slice and 400 for the high-2 slice. The Q- values are reported 



in units of /i ^ Mpc^. 



Value of 


Low-z 


slice 


High-z 


slice 


^max 


h 


Q 


h 





^orig 


1.97 ±0.05 


21.7 ±2.6 


1.98 ±0.03 


17.1 ±1.5 




2.03 ±0.07 


16 ±8 


1.96 ±0.04 


21 ±5 


iorig/4 


1.99 ±0.12 


33 ±45 


2.00 ±0.07 


-12 ±24 



the threshold luminosity varies slowly across that range, 
so we expect the bias to not vary significantly across the 
redshift range. This expectation has been confirmed in 
previous angular clustering studies which found ~ 15% 
variation from z = 0.2 to z = 0.6 [t^I, and also by our 
own bias analysis which finds no significant difference 
between the two bins. Thus we conclude that for the 
purposes of ISW work (where we have a ^ 1.3(2.7) sigma 
signal for low-z LRG (high-z LRG) correlation), variation 
of the LRG bias within an individual photo-z bin (0.2-0.4 
or 0.4-0.6) can be neglected. 



We calculate the possible contribution from magnifi- 
cation bias given the redshift distribution of the LRGs 
and also an assumed cosmology. We find that the pos- 
sible contribution from magnification bias is 100 — 1000 
times (depending on the scale) smaller than the actual 
signal. Therefore magnification bias is not contributing 
significantly to our signal. 



C. SDSS quasars 

The function fi{z) for the quasars is more uncertain 
than for the LRGs. This is in part due to the Hmited 
spectroscopic coverage available, but also the difficulty 
of constructing quasar photo-zs and the lower clustering 
amplitude, which leads to noisier estimates of bias pa- 
rameters. The basic procedure for obtaining /i(z) is thus 
to find a region of sky with as high spectroscopic com- 
pleteness as possible while still retaining a large area; use 
this to obtain a preHminary estimate n(z); and then fit 
for the bias parameters using clustering data, of which 
several are needed if n(z) is multimodal. The remain- 
der of this section describes the details of the /i(z) de- 
termination and what possible errors can be introduced 
by spectroscopic incompleteness, stellar contamination, 
redshift-dependent bias, and cosmic magnification. 

In order to determine the redshift probability distri- 
bution, we began by constructing a set of five rectangles 
that He within the cov erag e area of the SDSS, 2QZ [7ll |. 
6QZ [nl, and 2SLAQ [zal surveys. These rectangles He 
along the equator (the declination range is — 01°00'36" 
to ±00°35'24") and cover the five RA ranges 137-143°, 
150-168°, 185-193°, 197-214°, and 218-230°. There is 
a significant amount of area with coverage from all sur- 
veys that is rejected as it was found to have lower com- 
pleteness in 2SLAQ because there is less plate overlap. 
Spectra in SDSS were required to have high confidence 
(zConf > 0.95) [Hi and those in 2QZ, 6QZ, and 2SLAQ 
were required to be of high quality (quality== 11) [7l] |. 

Our coverage rectangles contained a total of 1410 low- 
redshift and 1269 high-redshift photo-z quasars; these 
numbers are lower than the product of the spectroscopic 
coverage area and the number density of photo-z quasars 
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because some parts of the latter catalogue were rejected 
by our stellar density cuts. Of the low-redshift photo-z 
quasars, we found that 257 (18%) had no spectroscopic 
redshift determination or low quality ones, 58 (4%) were 
identified as stars, and the remaining 1095 (78%) are ex- 
tragalactic. For the high-redshift sample these numbers 
are 208 (16%), 13 (1%), and 1048 (83%) respectively. 
From this data we construct a preliminary redshift proba- 
bility distribution npi.eiim(z) for each of the photo-z slices 
using a kernel density estimator, 

1 1 

np.cii„.(z) - E^^^"''"'''''^'''' (25) 
-'Vex \/ Ztt a 

where TVox is the number of matches to extragalactic ob- 
jects, Zk is the redshift of the fcth object, and a is the ker- 
nel width. The estimator is consistent in the limit that 
the number of objects TVex —>■ oo and cr —> at fixed A^cxC. 
In practice, a must be chosen to be small compared to 
the width of any real features in the redshift distribution 
(otherwise these are artificially smoothed out) , and large 
enough to smooth out shot noise (and redshift-clustering 
noise, if significant). We have used a = 0.04 (using 
a = 0.02 changes the fit bias by only 5%). This prelimi- 
nary distribution is shown in the top panel of Fig.[9l The 
redshift distributions in the two photo-z quasar slices are 
multimodal due to the nature of the photo-z error dis- 
tribution: the quasar spectra redward of Lyman-a are 
usually characterized by a roughly power-law continuum 
with superposed emission lines. This means that quasar 
colors oscillate as emission lines redshift into and out 
of the SDSS filters, resulting in an (approximately) self- 
intersecting locus in color space and many degeneracies 
in the photo-z solution. 

If the quasar bias were constant and magnification 
bias negligible, then we would have simply /i(z) = 
cnproiim(^:), with the proportionality constant c being the 
product of the bias and the probability for a photo-z 
quasar to actually be extragalactic. This constant could 
then be determined by fitting the amplitude of the quasar 
autocorrelation function, as has been done in most past 
ISW studies. However, in the real Universe quasars are 
known to have an evolving bias, which is potentially sig- 
nificant across the redshift range considered, and at red- 
shifts z ^ 0(1) lensing magnification can become sig- 
nificant. The magnification can be calculated from the 
slope a of the quasar counts near the 5 = 21 magnitude 
limit, which gives a = 0.82 for the low-z sample and 
a = 0.90 for the high-z sample. In principle the cut on 
the M-band magnitude error ((t„ < 0.5) could have an ad- 
ditional effect since magnification will reduce (t„ ; however 
this is not an issue for us since at the g = 21 threshold, 
for UVX objects we will have m < 22 where the typical 
magnitude error is < 0.5 even accounting for extinction 
(^M,max = 0.26). Since for these samples a — 1 is small, 
we compute the magnification bias using npi.oiim(z) in 
place of the true distribution n(z). That is, we replace 
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FIG. 9: (a) The preliminary quasar redshift distribution, con- 
structed from the successful matches to spectroscopic data, 
(b) The best-fit f{z) for the two quasar samples as described 
in the text for the WMAP cosmology. 

Eq. {21]) with 

/,(z)«6,(z)n,(z) 

. prelim (z')dz'.(26) 

This leaves only the problem of constraining the prod- 
uct bi{z)Ili{z) using the clustering data, i.e. the quasar 
power spectrum and quasar-LRG cross-power. Unfor- 
tunately the data is not capable of constraining a full 
model-independent distribution, so instead we write 

b,{z)n,{z)D{z) = A{z)U,^ 

prelim 

(z), (27) 

where D{z) is the growth factor, and A{z) is a piecewise 
constant function of z. This is equivalent to assuming 
that the clustering amplitude (divided by spectroscopic 
completeness) of the quasars is constant in redshift slices, 
which has been found to be a better a ppr oximation than 
constant bias in most quasar surve ys |73 |. For compari- 
son, the empirical "Model 3" of Ref. Q predicts b(z)D{z) 
to change by only 5% from z = 0.65 to 1.45, and by 
13% from z = 1.45 to 2.00. For the more recent model, 
Eq. (15) of Ref. [zi, these numbers are 24% and 15% 
respectively. At higher redshifts (z > 3) there is a sharp 
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increase in b{z)D{z) [t^ but UVX-selected samples do 
not contain objects from this redshift range. 

We constrain A{z) in as many redshift sHces as can be 
constrained using the data. In particular since the quasar 
redshift distributions are multimodal, we would Hke to be 
able to fit a different clustering amplitude in each peak. 
The treatment of the two quasar samples is slightly differ- 
ent due to the availability of different information in their 
redshift ranges, so we now discuss their redshift distribu- 
tions separately. In each case, the autopower spectra 
were fit to linear theory up to Z = 160 {k = 0.1/iMpc~^ 
at z = 0.6) and the quasar-LRG cross-spectra were fit up 
to / = 140 (fc = 0.1/iMpc-i at 2 = 0.5). 



1. Low-z sample: 0.65 < Zphoto < 1.45 

For the low-z quasar sample, we can only constrain one 
redshift slice. An examination of Fig. [9] shows that the 
distribution is actually trimodal, with peaks at z = 0.32, 
1.24, and 2.20. A fit assuming a constant A yields A = 
1.36±0.10, with xVdof= 36.32/27 {p = 0.11). Almost all 
of the weight for this comes from the central (z = 1.24) 
peak. We also ran two-slice fits to determine whether the 
clustering data constrain the amplitudes of the low- and 
high-redshift peaks. The first such fit is of the form 

which allows the low-redshift slice to vary (z = 0.52 is the 
local minimum of pieiim)- This fit gives Ai = 4.74 ± 
2.12 and A2 = 1.35 ± 0.10, with x^/doi= 33.77/26. We 
also tried a two-parameter fit in which the high-redshift 
slice is allowed to vary: 

^(") = \A^ z>1.83 (29) 

(the local minimum of Ili^proiim between the main and 
high-redshift peaks is at z = 1.83). This fit gives 
A[ = 1.37^^9 and A'^ = 0.0 ± 8.7 (Icr), with xVdof= 
36.31/27. The errors on A[ are highly asymmetric in this 
case because the constraint comes mainly from the quasar 
autopower; A[ and A'2 are then degenerate because one 
only knows the total power, not how much comes from 
each redshift sHce. The shape of the power spectrum 
breaks this degeneracy in principle, however in practice 
it is far too noisy. The fact that the high-redshift slice 
cannot give negative power accounts for the "hard" upper 
limit on Aj. 

From this exercise we conclude that the clustering data 
cannot independently measure the bias in either the low- 
or high-redshift peak. The reasons are different in each 
case. The low-redshift peak contained only 1.7% of the 
spectroscopic identifications, and thus almost certainly 
contains only a very small fraction of our quasars. This 
peak lies at the same redshift as the low-z SDSS LRGs, 



and the quasar-LRG cross-correlation is the major con- 
straint on Ai. Unfortunately this cross-correlation is 
drowned out by the enormous Poisson noise contributed 
by the quasars in the other two peaks, and is detected 
at only 2.2fT. On the other hand, the LRGs oversample 
the cosmic density field on linear scales and cover the 
same region of sky as the quasars. One would thus ex- 
pect that since the LRG-quasar correlation is only seen at 
this low significance, and the ISW effect from this redshift 
range contributes only a small fraction of the power in 
the CMB, the contribution of the low-redshift peak to the 
quasar-ISW correlation would be statistically insignifi- 
cant. We find that the predicted peak of the quasar-ISW 
l{l + l)Ci/2-K for only the low-redshift peak quasars is 
lower than the entire sample (high-z QSO) by O.OlS^if, 
which is significantly smaller than the error on the cross 
correlation. This is run using a WMAP-3yr parameters. 

The high-redshift peak contains 10% of the quasars. 
Its amplitude must be measured in autocorrelation due 
to the lack of other samples at that redshift, which is 
a serious drawback since only 1% of quasar pairs come 
from the high-redshift peak. An alternative approach 
to constraining its amplitude would be cross-correlation 
against the spectroscopic quasar sample at 2.0 < z < 2.5, 
but we did not pursue this approach here. 



2. High-z sample: 1.45 < Zphoto < 2.00 

The high-z photometric quasar sample also has a tri- 
modal distribution: there is one peak at z = 0.22, a 
second ai z ~ 0.58, and a third at z = 1.80. In this 
case however, it is the highest-redshift peak that contains 
most of the objects, with the middle peak in second place 
and only a few objects in the lowest-redshift peak. This 
situation makes it both possible and necessary to fit sep- 
arate amplitudes for the peaks; in this case we will find 
that two amplitudes can be constrained, one for the two 
low-redshift peaks and one for the main (high-redshift) 
peak. 

As a first step, we attempt to fit all three of the peaks 
with separate amplitudes, 

( Ai z < 0.33 
A{z) = I A2 0.33 < z < 1.18 . (30) 
[ A3 z>1.18 

This leads to the results Ai = 8.2 ± 4.5, A2 = l.Mtl^l, 
and A3. = 1.38tnt {la), with xVdof= 23.58/25. The 
large error bar on Ai indicates that this parameter can- 
not be constrained from the data, so we instead try a two- 
slice fit in which we fix Ai = A^. This fit gives the tighter 
constraints ^1=^2 = 1.59 ±0.61 and A3 = 1.35 ±0.10, 
with x^/dof= 25.75/26, and it is what we use for the rest 
of the paper. 
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3. Redshift Distribution Summary 

The quasar autopower spectra and quasar-LRG cross- 
spectra, along with the model fits, are shown in Fig. [TOl 
For the QSOO sample, there is excess power 3cr above 
the prediction) in the lowest-^ bin, corresponding to a 

2% RMS fluctuation in the number density on scales 
of ^ 30 degrees. The two most obvious sources of such 
power are stellar contamination and photometric cali- 
bration errors. Given that ~ 5% of the photometric 
"quasars" are actually stars [HHl and that the relative 
photometric calibration across the sky in SDSS is esti- 
mated to be ~ 2% in the u band (the worst band, but 
one very important for quasar work) [tOI, either of these 
seems plausible. In any case, these very low multipoles 
were not used in fitting the redshift distribution in either 
auto- or cross-power. 

It is essential to test the robustness of the quasar fits, 
in particular against the possibility of nonlinear cluster- 
ing affecting the range of multipoles used in the fits. The 
first way we do this is by repeating our analysis using the 
nonlinear matter power spectrum of Smith et al. [62| in 
place of the linear power spectrum. In the analysis with 
the nonlinear spectrum, the amplitude A for the low-z 
quasar slice increases by +0.02, and the amplitudes for 
the z < 1.18 and z > 1.18 parts of the high-z quasar slice 
increase by -1-0.08 and -1-0.02, respectively. If we restrict 
our attention to the lowest multipoles / < 100 (instead of 
cutting at 140 or 160), these changes are -1-0.02, —0.14, 
and -1-0.03. In each case the change is very small com- 
pared with the error bars. Thus we do not believe that 
nonlinear clustering is affecting our fQso{z) estimates. 



D. NVSS 

The function f{z) for NVSS is the hardest to obtain 
because there are no spectroscopic samples of NVSS ob- 
jects that have sufficiently high completeness to obtain 
the redshift distribution. Past ISW analyses with 
the NVSS have been based on the radio luminosity func- 
tion $(i, z) of Dunlop & Peacock [t^I, which itself was fit 
to a combination of source counts, redshifts for some of 
the brightest sources, and the local luminosity function. 
A constant bias was then assumed. The redshift distri- 
bution so obtained is reasonable, however it has three 
major drawbacks: (i) the redshift probability distribu- 
tion n(z) for the faint sources (which make up most of 
the sample) is constrained only by the functional form 
used for the luminosity function and not by the data; 
(ii) it does not give the redshift dependence of the bias, 
which could be very important since the redshift range 
is broad, and the typical luminosity of the sources varies 
with redshift; and (iii) the absolute bias b is constrained 
using the NVSS autopower spectrum, which is known to 
contain power of instrumental origin and hence is prob- 
ably a less reliable constraint than the cross-correlation 
against other surveys. The alternative method to mea- 



TABLE IV: Details of the cross-correlation of NVSS with the 
eight other samples. The second and third columns show 
the fraction of objects in each of the samples that match to 
the NVSS, i.e. rii,Nvss/'T-i. Results are presented for two 
matching radii, 40 and 20 arcsec. The final two columns show 
the range of multipoles used in the cross-correlation. 



Sample 


rii.Nvss/iii 


Multipoles used 




40" 


20" 


^min 


^max 


2MASS 12.0 < 7^20 < 12.5 


0.1317 


0.1302 


10 


14 


2MASS 12.5 < 7^20 < 13.0 


0.0802 


0.0787 


10 


14 


2MASS 13.0 < K20 < 13.5 


0.0473 


0.0455 


10 


24 


2MASS 13.5 < K20 < 14.0 


0.0292 


0.0280 


10 


36 


SDSS LRG \ow-z 


0.0450 


0.0425 


10 


87 


SDSS LRG high-2 


0.0263 


0.0249 


10 


139 


SDSS QSO \ow-z 


0.0180 


0.0192 


10 


239 


SDSS QSO high-z 


0.0189 


0.0207 


10 


159 



sure f{z) is by cross-correlation against the other samples 
whose redshift distributions are known. This method is 
adopted here, since it does not have any of the aforemen- 
tioned problems. Its main drawback is that the other 
samples only probe the range out to z ^ 2.6, and little 
data is available to constrain f{z) above that. 



1 . Procedure 

In order to measure the effective redshift distribution of 
NVSS, we must first obtain the cross-correlation of NVSS 
with each of the eight other samples (the four 2MASS 
samples, and two samples each of LRGs and quasars). 
This is done by using the same angular cross-spectrum 
estimation method as was used for the ISW analysis, and 
the cross-spectra are shown in Fig.[TTJ The main subtlety 
that arises is that the cross-spectrum C^"* (where i and 
j are LSS samples) can actually contain Poisson noise if 
there are objects that are in both samples. The Poisson 
noise term is of the form 



c;^ =c;^(Lss) 



niUj 



(31) 



where is the number of sources per steradian in catalog 
i, and fiij is the number of sources per steradian that 
appear in both catalogs. In order to measure fiy we must 
match the NVSS to each of the other samples. Note that 
the positional errors in NVSS are typically several arc 
seconds, and consequently there will always be some false 
matches. Therefore we estimate the fraction of matches 
as 



"-i.NVSS _ 



match 



7T0i 



,'T'NVSS, 



(32) 



where A^match is the number of matches within some ra- 
dius 6'ijiax, and A^j is the number of sources in catalog i 
in the NVSS mask. This was estimated for radii Omax of 
40 and 20 arcsec, and the results are shown in Table [iVl 
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FIG. 10: The model fits to the power spectra of the quasars and their cross-correlation with the LRGs. The low- and high-z 
quasar slices are denoted "QSOO" and "QSOl" respectively, and a similar nomenclature is used for the LRGs. The model fits 
using linear theory are shown with the solid lines over the range of multipoles used in the fit. The dashed lines show the 
extension of the model across the remaining range of multipoles. Note that for the highest multipoles the linear theory is 
expected to break down. 



We next computed the cross-power spectra between 
NVSS and each of the other samples. These spectra 
(after subtraction of the Poisson term) are shown in 
Fig. [iTJ The redshift distribution was then fit to the 
cross-power spectra. In this fit the minimum multipole 
used is ^min = 10 (below which there is a large amount of 
spurious power in the NVSS map) and the highest-Z bin 
used was determined by the formula ^max = fcmax^yi.20 5 
where fcmax — 0.1/iMpc~^ is the smallest scale to be fit 
and Da, 20 is the distance corresponding to the 20th per- 
centile of the window function for that sample as de- 
fined in Appendix El We have fit /nvss(^) with a F- 
distribution, 

/nvss(^) = St] , ^cff^"e— . (33) 
I (a) 

This function has three free parameters, bcS, z^, and a. 
Of these the normalization bas may be viewed as an ef- 
fective bias in the sense that / /nvssI^) dz = bcs; in the 
absence of cosmic magnification this would be the bias 
averaged over the redshift distribution. The peak of the 
distribution is at 2^, and a controls the width of the dis- 
tribution. The parameter fit gives 6off = 1.98, = 0.79, 
and a = 1.18. 



2. High-redshift tail 

The above analysis of the NVSS distribution involved 
cross-correlations against several samples at < z < 2. 



(The QSOO sample has a small number of objects at 
2.0 < z < 2.6, however they have no significant impact 
on the fitting of the QSOO x NVSS cross-spectrum.) Thus 
it leaves open the issue of whether there is a tail of ob- 
jects at high redshift, z > 2. Since /(z) is a product of 
bias times redshift probability distribution, it need not 
be normalized - / /(z) dz can have any value - so there 
is no way to tell from the cross-correlation analysis alone 
whether a portion of the sample is missing. If we also 
use the NVSS autopower spectrum then in principle one 
can determine whether an additional source of angular 
fluctuations is necessary. However the angular cluster- 
ing at flxed angular scale I is much stronger at low than 
high redshift, and the NVSS autopower spectrum is of 
low signal-to-noise ratio and possibly contaminated by 
systematics, so we have not chosen this strategy. 

An alternative approach to the high-z tail is to directly 
match against optical/ NIR catalogs. One can then use 
the rriK — z relation or (if multiband imaging is avail- 
able) photometric redshifts. There are always some ra- 
dio sources without optical identiflcations, however this 
method enables one to set an upper limit to the number 
of NVSS sources that can be at high redshift. For our 
analysis, we have matched against the COSMOS fleld, 
which has a modest solid angle (2 deg^), multiband imag- 
ing allowing good photometric redshifts, and deep high- 
resolution coverage with the VLA. Area is required due 
to the low density of NVSS sources (40 deg~^), and high- 
resolution radio images are required to uniquely identify 
an NVSS source with an optical counterpart due to the 
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FIG. 11: The cross-spectra of N'VSS with the other samples. The solid lines show the linear theory predictions in the region 
used for the fits, and the dashed lines show the extension to higher or lower multipoles. Note that for the highest multipoles 
linear theory is not valid. 



large positional uncertainty in the NVSS 
faint sources) 15811. 



7 arcsec for 



The COSMOS field contains 87 NVSS sources that 
pass our cuts. We began by matching these to the 
VLA-COSMOS observations, which are much deeper and 
have typical positional uncertainties of ^ 0.2 arcsec [t^ . 
Of the NVSS sources, 79 have a match within 30 arc- 
sec (we take the nearest source in the event of multiple 
matches). The 79 VLA-COSMOS sources that match to 
NVSS are then matched to the optical catalog [t^; there 
are 64 successful matches within 1 arcsec. This repre- 
sents 74% of the original NVSS catalog. It is of course 
possible that there are some false matches. By adding 
up nirO'^ for each NVSS source, where h is the density 
of VLA-COSMOS sources and 6 is the distance to the 
nearest VLA-COSMOS source (or 30 arcsec if the NVSS 
source had no match), we estimate that there are ^ 5 
false NVSS/VLA-COSMOS matches. A similar argu- 
ment suggests that ^ 0.5 false matches of VLA-COSMOS 
to the COSMOS optical/NIR catalog. Thus we expect 
that 58.5 of the matches are correct, corresponding to 



67% of the initial NVSS catalog. 

We show the photometric redshift distribution of the 
matches (according to Mobasher et al. [sOl) in Fig. [121 
Our best-fit /nvss(-z) (with the F distribution) has 24% 
of the bias-weighted source distribution at z > 2 and 8% 
at 2: > 3; if the source bias increases with redshift, as 
usually found for optical quasars, this number would be 
lower. From Fig. [12] we see that only 2 out of 64 matches 
fall at z > 2, i.e. the high-redshift tail of the F distri- 
bution can only exist in reality if (i) most of the 26% of 
the sources with failed matches to COSMOS optical/NIR 
data are actually at z > 2, or (ii) the sources at z > 2 
have a large bias. Both (i) and (ii) are physically plausi- 
ble but we have no direct evidence for them. 

The conservative solution in this case is to consider 
two Hmiting cases for the redshift distribution of the 
sources at z > 2. One case, which gives the minimal 
lensing signal for all cosmologies, and the minimal (max- 
imal) ISW signal for ACDM (closed) cosmologies, is to 
set /nvss = at z > 2. In the opposite limiting case, we 
have assumed that all failed and incorrect NVSS matches, 
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FIG. 12: The redshift histogram of NVSS sources matched 
to COSMOS using the Mobasher et al. [13] photometric red- 
shifts. The dashed line is the fit three-parameter /nvss(-z), 
normalized to unity (i.e. the redshift distribution assuming 
constant bias and negligible effect from magnification). 



and all sources with Zphoto > 2 (i.e. a total of 35%) are 
at z > 2, and have four times the clustering amplitude 
measured for the optical quasars (QSOl sample), e.g. 
b{z) = 4 X 1.35/ D(z) (where D is the growth factor) for 
the fiducial cosmology; the shape of /nvss(-z) at z > 2 
was left unchanged from the F-distribution fit. In order 
to understand the change of ISW and CMB-lensing sig- 
nals due to changes of our assumption of the high-z end 
of the redshift distribution of NVSS, we look at two dif- 
ferent redshift distributions, one with nothing at z > 2 
(minimal model) and the other with a "maximal" num- 
ber of sources (assuming clustering strength 4 times of 
the optical quasars and all the failed optical IDs are at 
z > 2). We find that the signals for both ISW (average: 
7.8%) and CMB-Lensing change by less than 10%, there- 
fore, one won't expect the unidentified high-z tail of the 
NVSS sources be a problem in our analysis. 



3. Constraints, robustness, and alternatives 

While the fit parameters are formally determined by 
the it is useful to graphically display the constraints 
in order to show what parts of the distribution are con- 
strained by which data. This we have done in Fig. [131 
For each of the eight samples, we have plotted on the 
vertical axis the constant /nvss value that provides the 
best fit to cross-correlation with that sample and its la 
error bar. The horizontal position is determined by the 
following procedure. We show in Appendix [A] that the 
estimated constant /nvss is actually given by an integral 
over some window function, 

/>oo 

(/nvss)=/ W{z)hyss{z)dz, (34) 
Jo 
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FIG. 13: The constraints on the NVSS redshift distribution 
from the cross-correlations with the other eight samples. The 
horizontal error bars show the redshift window functions as 
described in the text. The dashed line shows the result of us- 
ing the redshift distribution based on the Dunlop & Peacock 
[77|| luminosity function assuming constant bias and neglect- 
ing magnification, as has been done in most ISW studies. 



where the window function yV(z) integrates to unity. The 
horizontal position of the data points in Fig. [13] is the 
median of the window function, i.e. the redshift z where 
/q W{z') dz' = 1/2. The error bars extend from the 20th 
to the 80th percentile of the window function. 

Finally we wish to compare the redshift distribution 
we have obtained to that used in previous ISW studies. 
The previous results were based on the radio luminos- 
ity function of Dunlop & Peacock [t^I- In each case, it 
appears that the authors used the luminosity function 
and fc-correction based on the spectral index to infer the 
redshift distribution, assumed constant bias and negligi- 
ble magnification, and determined the one free parame- 
ter (the bias) by fitting to the autopower spectrum. If 
we do this using the fiducial WMAP cosmology and our 
autopower spectrum we find b ~ 1.7, and the function 
/nvss (2) = OT(z) obtained is shown as the dashed fine 
in Fig. [131 This curve, while roughly consistent with the 
NVSS-quasar and NVSS-LRG correlations, badly over- 
predicts the NVSS-2MASS correlation. Note that the 
problem cannot be fixed by changing the single bias pa- 
rameter: if h were reduced by a factor of ~ 3 to fit the 
2MASS data, then the LRG and quasar data would be 
discrepant. 

There are several possible explanations for this: 

1. The shape of /nvss(-z) is being modified by magni- 
fication bias. 

2. The extrapolation of the luminosity function to 
faint sources at high redshift by Dunlop & Peacock 
is in error. 
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3. It is possible that the Dunlop & Peacock redshift 
distribution accurately describes the NVSS sources, 
but the bias increases with redshift so as to produce 
the shape seen in Fig. [131 

4. The cut imposed by us (and by other ISW groups) 
that requires NVSS sources to be unresolved is se- 
lecting against nearby objects, and hence pulling 
down the low-z part of the /nvssI-z) curve. 

Of these, possibility number [T] is easy to rule out. Ap- 
plication of Eq. ((2TI) impHes that /nvss(^) has a max- 
imum change due to magnification bias of 0.09|a — 1| 
{z = 0.55), and a smaller change at lower redshift 
(0.03 1 a - 1\ at z = 0.1), where a = -dlogN/dlogF 
is the source count slope. The NVSS point source counts 
suggest a slope of 0.99 between 2.5 and 5.0 mJy, and 0.95 
between 5 and 10 mJy, which suggests that the effect of 
magnification bias on A/nvss(-z) is at most of order 0.01. 
In order to accommodate the discrepancy of A/nvss(2;) 
between our result and the Dunlop & Peacock distribu- 
tion of ~ 0.6 at z < 0.1, we would need an absurd slope, 
a « -20. 

Distinguishing among the remaining three possibilities 
is harder. We believe possibility number [2] is unlikely 
because the discrepancy between Dunlop & Peacock and 
our work occurs at low redshift where their luminosity 
function should be most reliable: this regime is con- 
strained by the local source counts rather than by extrap- 
olation. Redshift-dependent bias (possibility number [3]) 
exists for most samples of objects and there is no rea- 
son to expect it to be absent for NVSS. However, based 
on the Dunlop & Peacock dN/dz and our /nvss(^); the 
bias would have to change from ~ 0.4 at z = 0.1 to ~ 2 
at z = 0.5. Such a large variation, combined with the 
unusually low value of the bias at z = 0.1, suggests that 
this is not the full explanation. The final possibility (|4]) 
is the removal of extended sources. This is hard to as- 
sess because of the low density of extended NVSS sources 
above our fiux cut (^ 8deg~^). Of the 20 such sources 
in the COSMOS field, 19 match to VLA-COSMOS and 
13 of these matches are found in the COSMOS opti- 
cal/NIR catalog. It is worth noting that 8 of these (62%) 
have Zphoto < 0.5, versus 30/64 (30%) for the unresolved 
NVSS sources. This appears to go in the right direction, 
however it is difficult to make quantitative statements 
about whether the extended sources actually resolve the 
discrepant redshift distributions because of the unknown 
(but probably large, especially for the low-z part of the 
distribution) sampling variance error bars. 

In summary, while the full explanation for the differ- 
ence between our /nvss(-2) and that of Dunlop & Pea- 
cock remains unknown, it seems Hkely (based on process 
of elimination) that a combination of redshift-dependent 
bias and our rejection of the unresolved NVSS sources 
plays a role. Magnification bias is ruled out as the ex- 
planation, and the discrepancy occurs in a regime where 
the extrapolations used in Dunlop & Peacock probably 
do not matter. 



VI. SYSTEMATICS 

We investigate various systematic effects in our corre- 
lations utilizing a specific multipole range. We choose 
these multipole bins based on two criteria. First, they 
should not be affected by non-linearities. Second, they 
should not be affected by any of the systematic effects 
in a significant way. We therefore only utilize the mul- 
tipoles corresponding to fc < 0.05/iMpc~^ and we also 
discard the first ^-bin for all samples since it is affected 
by the galactic foreground contamination. The specific 
Z-bins that are utiHzed are tabulated in Table IVl 



A. Dust Extinction 

Since it is possible that incorrect dust extinction sys- 
tematically adds signals to our ISW cross correlation, 
we cross correlate the reddening maps jl^l in the same 
manner as we cross correlate each of our sample to the 
cosmic microwave background. If there is a systematic 
effect contributed via dust extinction, it will show up as 
a correlation, we can then estimate the effect and correct 
it from our tracer-cmb correlation. 

In order to the verify that dust extinction does not 
affect our results, we constructed a vector f of the esti- 
mated spurious cross-spectra AC^^ . The spurious cross- 
spectra were computed by taking the cross-power spec- 
trum of the CMB with the reddening map and multiply- 
ing by an estimate of dSg/dE{B — V). Note that f has 
an entry for each £-hin for each sample, so it has a total 
length of 42. We then compute the quantity (the deriva- 
tion of this quantity and its relevance to understand con- 
tamination from extinction is detailed in Appendix ICj) : 

E,,t = f^C-if. (35) 

Here C is the total 42 x 42 covariance matrix that is gen- 
erated using looking at the covariances of the correlation 
with each tracer sample and the Monte-Carloed CMB 
temperature map (the MCI procedure in the terminol- 
ogy of Cabre et al. [sil; see Sec. lVIIB 2l for details). 

Here y/E^xt is the maximum number of sigmas at 
which the effects of dust extinction could be detected if 
we knew all cosmological and redshift distributions per- 
fectly; if Eext ^ 1 then the dust extinction cannot have 
any statistically significant effect on any quantity derived 
from the cross-power spectrum, including cosmological 
parameter estimates. We estimate that dSg/dE{B — V) 
= -0.1 (all 2MASS samples). For the SDSS samples we 
did a Poisson-weighted fit to the LRG and quasar over- 
densities versus E{B — V) (see Fig. [2]); this gives —0.76 
(low-z LRGs), -0.18 (high-z LRGs), -1.06 (low-z QSOs), 
and —0.26 (high-z QSOs). (The Poisson error bars are all 
within 2a of zero so there is no evidence that any of these 
derivatives is nonzero.) We ignore extinction for NVSS 
since it is at radio frequencies. This gives ^/E^xt =0.23, 
so the dust extinction is not having a significant effect. 
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B. Galactic foregrounds 

To test whether galactic foreground contamination is 
important in our analysis, we cross correlate the tem- 
plates of Galactic emission with the tracer overdensity 
maps. The galactic foregrounds that must be considered 
in producing a template at higher frequencies are free- 
free and thermal dust emission; at lower frequencies an 
additional component is present whose physical origin re- 
mains uncertain but which may include hard synchrotron 
emission [3| or spinning or magnetic dust [82, We 
have used Model 8 of Schlegel et al. [33|, Finkbeiner et al. 
[s^ l for thermal dust and the Ha line radiation template 
of Finkbeiner [s^l rescaled using the conversions of Ben- 
nett et al. [3| for free-free radiation (see js^l for further 
details). We then construct these maps in the same way 
as in WMAP temperature maps. Cross correlations be- 
tween these templates with each of the tracer overdensity 
maps are then performed. 

To understand the foreground contamination to our 
result we compute as above 

Efa = f/gC-if/g (36) 

where f/g is the vector of cross-power spectra of the LSS 
and foreground maps, and C is the Monte Carlo covari- 
ance matrix. Calculating the \jEfg we find that the low 
multipoles of some of the low redshift samples correlates 
with the galactic foreground. We investigate this further 
and realize that there is incidentally a low redshift clus- 
ter at low latitude, thus correlating with the foreground 
map. We therefore restrict our 1-range that contributes 
to our signal by leaving out the first multipole bins for 
all sample. For the remainder we get y^Efg = 0.66. 

C. Thermal SZ effect 

The thermal Sunyaev-Zeldovich (tSZ) effect has a rel- 
atively weak frequency dependence compared to the 
Galactic foregrounds, so we constrain it from theoreti- 
cal models. We look at the tSZ signal using the halo 
model, separating the effect of the tSZ signal into 1-halo 
term and 2-halo terms. 

The 1-halo term stands for the situation when the fiux 
added towards the CMB map via tSZ effect comes from 
the same halo as the one that hosts the galaxies that we 
are correlating them with. The theoretical prediction for 
the 1-halo term is: 

/AT P 
dF—n2DiN,F) (37) 

where TV is the number of galaxies in that halo, F is the 
fiux from the halo, ifg is the average number of galaxies, 
n2D{N,F) is the number of halos with N galaxies and 
fiux between F and F + dF . We then turn Eq. (|37|) into 
integrals over halo mass and comoving distance: 

Cf^lh)^ I % [ dM^0(M)^F(M,x), (38) 
J J Po n-g 



where 4'{M) is the fraction of the mass in haloes between 
M and M + dM, N{M) is the mean number of galaxies 
in a halo of mass Af , and F is the fiux from a halo of 
mass M at comoving distance x- 

The 2-halo term stands for situation when the fiux 
(from tSZ) comes from a different halo which hosts galax- 
ies that cross-correlate with the fiux. It is 

//V F 
dF—n2D{N)n2D{F)Ct{N; F), 

(39) 

where n2D{N) is the number of halos with N galaxies 
per steradian, n2D{F) is the number of halos with fiux 
between F and F + dF per steradian and Ci{N;F) is 
the cross-power spectrum between halos with TV galax- 
ies and those with fiux F. We then turn the Eq. l(39|) 
into integrals over the mass functions and cosmological 
distances: 

Cf^{2h) = [ % [ dM—(j>{M)b{M) 
J J Po 

xfix)PHnik)FiM,x), (40) 

where Piin{k) is the 3-D linear matter power spectrum. 

Now, what is left for us to do is to figure out what the 
fiux F is for tSZ effect. One should note that this method 
is not limited to the tSZ effect prediction, but any kind of 
correlations between galaxy number overdensity and fiux 
of any kind associated with the halos. For tSZ effect, the 
flux is 

^ = 2f'rcMB^^^^, (41) 

where r' is the mean Thomson optical depth per unit 
comoving distance, Tcmb is the observed averaged CMB 
temperature, /icm is the baryon fraction in the intraclus- 
ter medium, ft is the cosmic baryon fraction, fcs is the 
Boltzmann coefficient, Te{M) is the average temperature 
of electrons inside halos of mass M, rrie is the mass of 
electrons, c is the speed of light. 

In order to assess the effect of tSZ on the ISW corre- 
lation, we calculate the Cf^{lh) and Cf^{2h) with a 
high (Tg (0.92) in order to give a conservative estimate. 
We must also estimate N{M). For the 2MASS samples, 
we use N{M) of the satellites and the conditional lumi- 
nosity function from [s^l while assuming that there is 1 
BCG per cluster. This is a conservative estimate as some 
of the BCGs may fall out of the flux limit. For the LRGs, 
we use N{M) from [s^ for our calculation without mod- 
ification, as we use the same galaxy sample. The quasars 
and NVSS are both examples of active galactic nuclei, 
and are generally found in haloes of some mass range 
with a small probability [i.e. {N){M) < 1] usually inter- 
preted as the duty cycle. For these cases, we first obtain 
the redshift distribution (dN/dz) and bias. For NVSS, 
we assume that bias oc 1/D{a) where D{a) is the growth 
factor of scale factor a, as there is no better available 
information (our determination of /(z) is not capable of 
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separately distinguishing the bias from the redshift distri- 
bution). From the bias, we constrain the minimum halo 
mass that will host a QSO or NVSS object, and then ob- 
tain the duty cycle based on dN/dz. Duty cycles cannot 
exceed unity, so we cap fduty at 1 and above this use 
dN/dz to get minimum halo mass. Then, N{M) = fduty 
if M > Mmin and otherwise. 

We assess the level of contaminations by calculating 



Etsz = Cf^{lh + 2h)CCf'^{lh + 2h), (42) 



\ntsz I 



which is the tSZ analogue to Eq. 



35]). We find that 
0.109 using the i'-bins that are tabulated in 



Table |V] and thus thermal SZ effect is not a significant 
contamination for the ISW effect. 

We present our results for the tSZ contamination for 
the /-bins that we use in our analysis the cosmological 
parameter estimation in Table IVl 



D. Point source contamination 

Point source contamination is one of the main concerns 
that we have for analysis for cross correlation of CMB 
with large scale structure, as point sources add to the 
CMB, while they are probably correlated with the tracers 
of large scale matter density field. Therefore, we estimate 
the contamination from the point sources by estimating 
CYiy) by looking at the differences of cross correlation 
of the tracer samples with different frequency maps of 
WMAP. We estimate Cf at 61 GHz (V band): 



Ce{Ka)-CiiV) 



TKaV 



Ka 



ryv^, 



(43) 



where rx is the ratio of thermodynamic temperature to 
the antenna temperature of band X and we assume that 
T{v) is proportional to . We assess the level of con- 
taminations by calculating (similarly as above mentioned 
foreground analysis): 



TABLE V: The tSZ and point source contamination for each 
of the samples we used in the analysis. For tSZ the 1 halo 
and 2 halo terms are shown separately and combined. 



Sample 


I 


tSZ Ih 


[l{l + l)/2-K]Cf (mK) 
tSZ 2h tSZ l+2h 


pt src 


2MASS0 


6 


-0.0085 


—0.0458 


—0.0543 


-0.4056 


2MASS1 


6 


-0.0048 


—0.0324 


—0.0372 


-0.0743 


2MASS1 


11 


-0.0151 


—0.0574 


—0.0725 


0.0070 


2MASS2 


6 


-0.0027 


—0.0241 


—0 0268 


-0.0875 


2MASS2 


11 


-0.0086 


—0.0458 


—0.0544 


0.0216 


2MASS3 


6 


-0.0016 


—0.0182 


—0.0198 


-0.1717 


2MASS3 


11 


-0.0050 


—0.0375 


-0.0425 


0.0089 


LRGO 


18 


-0.0045 


—0.0196 


—0.0241 


0.0020 


LRGO 


31 


-0.0132 


—0.0394 


—0.0526 


0.0261 


LRGO 


43 


-0.0251 


—0.0574 


—0.0826 


0.0123 


LRGl 


18 


-0.0017 


—0.0064 


-0.0081 


0.0018 


LRGl 


31 


-0.0049 


—0.0173 


—0.0222 


-0.0379 


LRGl 


43 


-0.0094 


—0 0269 


—0 0363 


0.0109 


LRGl 


56 


-0.0159 


—0.0361 


—0.0520 


-0.0028 


LRGl 


68 


-0.0240 


—0.0460 


—0.0700 


-0.0332 


QSOO 


18 


-0.0003 


—0.0012 


—0.0015 


-0.0039 


QSOO 


31 


-0.0010 


—0.0036 


-0.0046 


0.0058 


QSOO 


43 


-0.0018 


—0.0067 


—0 0085 


-0.0254 


QSOO 


56 


-0.0031 


—0.0102 


—0.0133 


0.0097 


QSOO 


68 


-0.0047 


—0.0135 


—0.0182 


-0.0509 


QSOO 


81 


-0.0064 


—0.0164 


—0.0228 


0.0660 


QSOO 


94 


-0.0086 


—0.0193 


—0.0279 


0.0169 


QSOO 


110 


-0.0118 


—0.0230 


—0.0347 


0.0626 


QSOO 


130 


-0.0164 


—0.0278 


—0.0442 


0.1854 


QSOl 


18 


-0.0006 


—0.0010 


—0.0017 


0.0000 


QSOl 


31 


-0.0018 


—0.0027 


-0.0045 


-0.0169 


QSOl 


43 


-0.0035 


—0.0046 


—0.0081 


-0.0131 


QSOl 


56 


-0.0058 


—0 0068 


—0.0126 


0.0030 


QSOl 


68 


-0.0088 


—0.0091 


—0.0179 


-0.0073 


QSOl 


81 


-0.0121 


—0.0112 


—0 0233 


0.0332 


QSOl 


94 


-0.0163 


—0.0134 


—0.0297 


0.0627 


QSOl 


110 


-0.0223 


-0.0158 


-0.0381 


0.0801 


QSOl 


130 


-0.0311 


-0.0184 


-0.0494 


0.0794 


QSOl 


150 


-0.0413 


-0.0207 


-0.0620 


0.0924 


QSOl 


170 


-0.0530 


-0.0232 


-0.0763 


0.0223 


NVSS 


6 


-0.0001 


-0.0007 


-0.0008 


-0.0398 


NVSS 


11 


-0.0003 


-0.0020 


-0.0023 


-0.0124 


NVSS 


20 


-0.0010 


-0.0050 


-0.0059 


-0.0111 


NVSS 


31 


-0.0023 


-0.0091 


-0.0113 


0.0014 


NVSS 


43 


-0.0043 


-0.0135 


-0.0178 


0.0103 


NVSS 


56 


-0.0073 


-0.0179 


-0.0252 


0.0025 


NVSS 


68 


-0.0107 


-0.0217 


-0.0324 


-0.0141 



E^, = Cf (i^)C-iCf (i/) (44) 

We find that yjEps = 0.495 using the ^-bins that are 
tabulated in Table [Vl and thus point sources is not a sig- 
nificant contamination for the ISW effect. Note that this 
includes some effect from Galactic foregrounds (which 
probably dominate the low £'s), since any foreground ef- 
fects that have frequency dependence will show up in 
Cf'^{V). In particular one would be double-counting the 
Galactic foreground if one added Eps and Efg. 

We present the point sources contamination for £-bins 
we use for our analysis in the last column of Table [Vl 



VII. COSMOLOGICAL IMPLICATIONS 

A. Significance of ISW detection after rejecting 
contaminating bins 

After investigating all the listed systematics and tak- 
ing into account of the non-linearities, we decide to only 
take the £-bins as are listed in Table El The high-^ 
bins are cut off due to the non-linearities; we cut off 
all the bins that at the median redshift for the fidu- 
cial cosmology correspond to fc > 0.05/i Mpc~^ using 
k — + l/2)/r. This is a more conservative cut than 
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the usual k = 0.1/iMpc~^ but it must be remembered 
that in linear theory the ISW effect is sensitive to the 
derivative of D{a)/a which contains a cancellation from 
the growth of structure in the numerator and the scale 
factor in the denominator. Therefore nonlinear effects 
could be larger than one naively expects. We cut off the 
first ^-bin for all samples as these are most affected by 
Galactic foregrounds. 

We calculate the significance of each of the sample by 
the standard method. First, we compute the amplitude 
of the signal (Appendix El in our case, fiducial model is 
based on the WMAP 3- year parameters): 



A 



data 



\ ^theory 



C theory f-^_\^theory ' 



(45) 



where C*'^'^°'^^ is the vector of predicted cross-power spec- 
tra for the fiducial cosmology, Cf"'*"' is the vector of ob- 
served cross-spectra, and is the inverse-covariance 
matrix. We obtain by Monte Carlo simulation as 
described in the next section. 

The error is similarly computed with: 



1 



^ ^theory ^^-i^theory 



(46) 



and the significance in sigmas is obtained by the usual 
calculation. A/ a. The result is shown in Table [Vll 

In Fig. [14] we plot the amplitude (A) and its error using 
covariance matrices and fisher matrices from the corre- 
lation of the tracer sample with WMAP V-band, com- 
puted with angular and redshift weighting optimized for 
WMAP3 model, together with theoretical predictions for 
three cosmological models (open, closed and fiat) to il- 
lustrate the constraining power on fix from ISW effect. 
Flat model is WMAP3 model and by definition its the- 
oretical prediction is ^ = 1 (see Eq. |45]). The other 
two models were chosen to He along the WMAP degen- 
eracy curve (which essentially keeps fixed flmh^, ^bh? 
and 6*, defined to be 100 times the ratio of the sound 
horizon to the angular diameter distance to recombina- 
tion), although this does not imply they are necessarily 
good fits to the WMAP data: the ISW signal in the 
CMB power spectrum itself can break the degeneracy 
between the parameters that keep the angular diameter 
fixed, but because ISW is a subdominant contribution to 
primary CMB even on the largest scales its power to dis- 
criminate among models is Hmited. We can see that the 
predicted amplitude of ISW signal for ACDM is positive 
(using the standard sign convention) because at late time 
when cosmological constant becomes important growth 
of structure is decreasing in time relative to Einstein- 
de Sitter (EdS) model and the associated gravitational 
potential, constant at high redshift when the Universe is 
effectively EdS, begins to decay. The decay is larger if we 
decrease (for which we need to go to a sHghtly open 
universe to preserve angular diameter distance) , which in 
turn increases ISW. On the other hand, a closed universe 



with Q.m > 1 accelerates the growth of structure relative 
to EdS, so potential is growing and this model predicts 
ISW signal with opposite sign. While the sign is essen- 
tially determined by the growth rate, its amplitude and 
scale dependence depend on other cosmological parame- 
ters as well and vary as a function of redshift, as shown 
in Fig. [H 

As can be seen from Fig. [14] and Table IVII we have a 
detection of ISW signal in a number of data sets. Most 
convincing are SDSS LRGl and NVSS, both at about 
3cr, followed by LRGO, QSOl and 2MASS3 at 1.2-1.5cr 
evidence. Remaining data sets have significance below 
Icr, although only one among them has negative signal, 
opposite to ACDM model predictions. The overall sig- 
nificance of detection with ACDM weighting is 3.7cr. We 
emphasize that while we use optimal weighting of data 
to maximize the signal by downweighting the scales and 
redshifts where we do not expect the signal, this depends 
somewhat on the assumed model, so the significance of 
detection can be somewhat affected by this. For exam- 
ple, we could instead of ACDM have used a model that 
predicts an upward feature at ^ = 30 that only occurs at 
redshift around 0.5, therefore taking advantage of the 3 
sigma excess power seen in LRGl at that scale (Fig [4]). 
Using this model would give high weight to that feature 
and would lead to a higher significance of the overall de- 
tection. Of course such aposteriori procedure is not really 
waranted, but it does highlight the difficulty of compar- 
ing the significance of detection among different analyses, 
which may have used different priors. This problem is 
exacerbated if cross-correlation function analysis is used, 
as in most of the previous work, because in that case a 
narrow feature in Fourier space would spread out to a 
broader feature in correlation function. 

While we find a 3.7a detection we also note that the 
observed ISW signal exceeds the predictions of WMAP3 
ACDM model by about 2(t, since the fit gives A — 
2.23 ± 0.60 relative to model prediction A = 1. The 
discrepancy is reduced if we change cosmological param- 
eters somewhat and this is explored further in the next 
subsection using MCMC analysis. 

To show that our results are consistent throughout dif- 
ferent bands in WMAP, thus there is no significant con- 
tamination from frequency dependent systematics, we 
show the amplitude of ISW signal and associate one 
sigma error relative to the WMAP3 model for each of 
the sample for all of the WMAP bands (except K band) 
in Table (VUl The differences in frequency A(Q)-A(V) 
and ^(W)-yl(V) are all < 0.25cr and most are < 0.15cr, 
and there is no consistent sign. This reassures us that the 
frequency-dependent foregrounds are subdominant to the 
statistical errors in these higher-frequency bands. The 
comparison with Ka band, i.e. A{Ka)-A{Y), is worse 
especially for 2MASS0 (the difference is < 0.5cr for the 
other samples), probably due to Galactic emission. 
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FIG. 14: The ISW amplitude (A) and errorbars a{A) for all 
samples plotted along the redshifts compared with predictions 
of WMAP-3 year parameters. The fitting and errors for this 
figure used the Fisher matrices from the correlation of the 
tracer samples with the various WMAP maps. We also show 
the expected amplitude for 3 model Universes along the an- 
gular diameter distance degeneracy curve. We calculate the 
expected amplitude by substituting our observed correlations 
with predicted correlations for each of the model Universe 
and proceed in the same manner as described in Eg 1451 The 
three model Universes are: ACDM model with the WMAP-3 
year parameters (open triangle with dotted line); closed Uni- 
verse (open pentagons with short dashed line) 116=0.215, fl,n 
= 1.25, Qk = -0.29, h = 0.32, as = 0.61; and Open Universe 
(open squares with long dashed line) f2i,=0.015, fim = 0.089, 
Qk = 0.03, h = 1.20, as = 0.73. Note that the redshift distri- 
bution is very broad for NVSS, giving rise to the jump in the 
open model prediction, even though the effective redshift of 
the sample is nearly the same as for low redshift QSO sample. 



TABLE VI: Amplitude of ISW signal and the associate one 
sigma error relative to WMAP3 model and significance of 
detection for each of the sample and when we combine all 
samples. These are calculated using the covariance matrix 
that are derived from the correlations with the Monte Carlo 
CMB maps (as described in Ea l47I The overall signal is 2 sigma 
higher than WMAP3 model prediction. 



Sample 


Amplitude {A ± a) 


# sigmas 


2MASS0 


-2.01 ± 11.41 


-0.18 


2MASS1 


-f3.44 ± 4.47 


0.77 


2MASS2 


-f2.86±2.87 


l.GO 


2MASS3 


-f2.44 ± 1.73 


1.41 


LRGO 


+ 1.82 ± 1.46 


1.25 


LRGl 


+2.79 ± 1.14 


2.46 


QSOO 


+0.26 + 1.69 


0.16 


QSOl 


+2.59 + 1.87 


1.38 


NVSS 


+2.92 + 1.02 


2.86 


All Samples 


+2.23 + 0.60 


3.69 



B. MCMC methodology and Likelihood function 

1. MCMC methodology 

A major goal of this paper is to provide a full like- 
lihood function with which cosmological models can be 
compared to each other. Here we describe the details 
of the Hkelihood function construction and apply it to 
some simple cosmological model parametrizations. Our 
goal is not to give an exhaustive parameter estimation 
analysis, but just to provide some characteristic exam- 
ples of possible appHcations. We include both ISW anal- 
ysis of this paper and the lensing analysis of Paper IL 
However, the latter effect has small statistical signifi- 
cance and does not contribute significantly to the like- 
lihood analysis. We decided to test the following cos- 
mological models: flat ACDM model {VLmh? , ^bh'^, 0, 
T, Us, As), ACDM + flK (not assuming flatness), flat 
ACDM + w (assuming flatness, but allowing dark en- 
ergy to evolve) . Here fim is the matter density, fib is the 
baryon density in units of critical density, fix = ^K/H'^ 
is the curvature K expressed in terms of critical den- 
sity, h = _ffo/100km/s/Mpc is the Hubble parameter, 9 
is 100 times the ratio of sound horizon to angular di- 
ameter distance at recombination, r reionization optical 
depth and Ug and As are the slope and amplitude (at 
k = 0.05/Mpc) of primordial power spectrum. We also 
refit for the bias with the redshift distributions for each 
of the dataset used for each of the cosmological parame- 
ter sets which we calculate the for. There is a detailed 
description of the determination of bias and redshift dis- 
tribution in Section |Vl We Hmit our search to models 
with scalar fiuctuations only with no running of spectral 
index, no tensors, and no neutrino masses. We assume 
fiat priors on all of the parameters defined above. The 
priors we use are shown at Table [VTTTl In addition we im- 
pose 40km/s/Mpc < Hq < lOOkm/s/Mpc and that age 
of the Universe has to be at least 10 Gyr and at most 20 
Gyr. These priors are appHed to all the chains that we 
show in the paper (including those with WMAP alone). 

In most cases the intervals are sufficiently broad that 
the boundaries do not matter, with exception of WMAP 
only case with curvature or dark energy, where we apply 
additional prior with Hq > 40km/s/Mpc. We search the 
parameter space using COSMOMC 183 with likelihood 
function from WMAP 3 year analysis [6a] • We discuss the 
Integrated Sachs Wolfe likelihood function in the follow- 
ing section, and leave the discussion of the Weak Lensing 
of CMB likelihood function to Paper II. We test the con- 
vergence of our Markov chains following Dunkley et al. 



S. Integrated Sachs Wolfe likelihood function 

This section describes the ISW Hkelihood function. We 
utilize the amplitude from the galaxy-temperature cross- 
spectrum C^^ from cross correlating the CMB sky (V- 
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TABLE VII: Amplitude of ISW signal and the associated la error relative to WMAP3 model for each of the sample for the 
WMAP bands (i.e. Ka, Q, V, W). The fitting and errors for this table used the Fisher matrices from the correlation of the 
tracer samples with the various WMAP maps. 



Sample 




Amplitude A 






Ka 


Q 


V 


W 


2MASS0 


-9.04 ±8.21 


-3.54 ±8.19 


-2.01 ±8.11 


-3.38 ± 7.79 


2MASS1 


1.80 ±3.97 


2.73 ± 3.94 


2.17 ±3.93 


1.64 ± 3.86 


2MASS2 


2.16 ±2.66 


2.95 ± 2.65 


2.42 ±2.63 


2.04 ±2.61 


2MASS3 


1.74 ± 1.72 


2.56 ± 1.72 


2.58 ± 1.72 


2.39 ± 1.69 


LRGO 


2.00 ± 1.44 


2.05 ± 1.44 


1.86 ± 1.45 


1.92 ± 1.46 


LRGl 


2.67 ± 1.04 


2.59 ± 1.04 


2.85 ± 1.05 


2.92 ± 1.06 


QSOO 


0.62 ± 1.90 


0.39 ± 1.92 


0.61 ± 1.89 


0.63 ± 1.94 


QSOl 


2.41 ± 1.90 


2.17 ± 1.92 


2.36 ± 1.90 


1.93 ± 1.90 


NVSS 


2.56 ± 1.01 


2.80 ± 1.01 


3.04 ± 1.02 


2.88 ± 1.02 



TABLE VIII: The priors applied to the 3 different chains. 
Note that all priors are flat. 

Parameter minimum maximum 



for all models, 6 parameters 





0.005 0.1 




0.01 0.99 


e 


0.5 10 


T 


0.01 0.8 


Us 


0.5 1.5 


loge{10^°As) 


2.7 4.0 


for ACDM ± only 




-0.3 0.3 


for ACDM ± w only 


w 


-2.1 -0.1 



band) with the following samples: 2MASS (0-3), SDSS- 
LRG (low-z and high-z), SDSS-QSO (low-z and high-z), 
NVSS. When we construct the likelihood function, we 
need three items: (i) the "data", which is C^^ for each of 
the sample for each £-bin (ii) the theoretical prediction; 
and (iii) the covariance matrix of the {C|^}. 

The data vector consists of the measured C^^ in each 
£-h'm and for each LSS sample used. After our cuts there 
are 42 such bins remaining, when combining all samples, 
thus the data vector has length 42. 

We calculate this covariance matrix by first generating 
1000 simulated CMB skies of WMAP resolution and then 
cross-correlate each of the samples with these simulated 
CMB sky. We call these cf^""''^ . We then calculate 
the covariance among the samples by first calculating the 
{Cf^"""'^) by averaging over all the correlations with all 
the simulated maps, then we find: 

x(Cf'^-(Cf— ))). (47) 

Note that this is a 42 x 42 covariance matrix, and that it 
is not block-diagonal in the LSS samples because there is 
some overlap in sample volume. The Monte Carlo proce- 
dure, by considering many realizations of the CMB but 
the actual realization of the galaxies, includes the implied 



correlations between different LSS samples. 

The issue of how to construct error bars on estimates 
of the galaxy-temperature cross-spectrum Cf^ , or its 
real-space equivalent Wgrid), has been a contentious is- 
sue ever since the first claimed ISW detections were an- 
nounced. The methods used have ranged from Gaussian 
error estimates based entirely on the theoretical galaxy 
and CMB spectra, to jack-knife methods that are based 
entirely on the data. Among the intermediate options are 
the Monte Carlo approach used here (MCI in the termi- 
nology of Cabre et al. [sH) in which the real galaxy field 
is cross-correlated against many random reaHzations of 
the CMB. 

If we knew the CMB and galaxy power spectra per- 
fectly from theory or observation, we would like to use 
analytic Gaussian error estimates for Cf'^ or do Monte 
Carlo simulations of random CMB and galaxy fields. 
Unfortunately, the galaxy maps, particularly QSOO and 
NVSS, are subject to spurious power at large angular 
scales for which we have no good theory, and for which 
we cannot measure the power spectrum accurately due 
to sampling variance. However we do know the theoreti- 
cal CMB power spectrum so we can implement MCI. It 
would also have been possible (but computationally ex- 
pensive) to implement a jack-knife; we chose not to do 
so because of concerns that at low multipoles the jack- 
knife regions would not be independent Il4ll although we 
note that the Cabre et al. simulations [8l| suggest that 
at least in some cases this is not a significant problem. 
The MCI method is however subject to two biases that 
could understimate the errors: a "correlation bias" due 
to neglect of the galaxy-temperature correlation when 
determining the error bars, and a "realization bias" due 
to the fact that only one realization of the galaxy field is 
used. These biases are discussed in Appendix [BJ where 
we find them to be negligible. 

We construct the likeHhood function as the following: 

= K(obs)-(.T")(p)][C-i]^.K(obs)-(x'')(p)], (48) 

where is simply C^"^; the index ^ encodes both the 
i-hin and the sample used. We denote by a;^(obs) the 
observed correlations C^"^, and (a;'')(p) denotes the mean 
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FIG. 15: The predicted ISW signal for the low-z LRGs (above) 
and high-z QSOs (below) sample for sample open, closed, and 
flat ACDM models. Parameters are the same as in Fig. 1141 



value predicted for cosmological parameters p. Note that 
the vector is of length 42 and that all LSS samples 
are included in a single x^; we do not add the val- 
ues of different samples separately since they are corre- 
lated and such an addition would be invalid. Among 
the three components of the likelihood function, only the 
predicted C^^ needs to be re-calculated for each cosmo- 
logical model. 



C. Parameter fits 

We investigate the following cosmological models: (i) 
ACDM model {nth'^, n^h'^, 0, r, n,. As); (ii) ACDM 



model + Qk; and (iii) ACDM model + w. Note that 9 
is the ratio of the sound horizon to the angular diame- 
ter distance, while Ag is the the primordial superhorizon 
power in the curvature perturbation on 0.05/Mpc scale. 
The numerical results are shown in Table [Xl for both the 
full likeHhood (CMB+ISW+WL) and CMB alone. We 
also looked at the effect of WL (or ISW) separately in 
constraining cosmological parameters by analyzing a cos- 
mological model (ACDM + ^k) using only CMB+ISW 
(without lensing) . We find the constraints to be similar to 
the full case (CMB+ISW+WL), but with slightly larger 
errorbars (see table |IX| ). Note that for the CMB-only 
model including ^k, the Markov chain ran up against 
the Hq > 40 km/s/Mpc boundary, thus artificially tight- 
ening the constraints; this did not occur for the full 
CMB+ISW+WL chains. 

For the ACDM model, the combined constraints from 
WMAP+ISW+WL is only slightly improved over using 
WMAP alone, but does lead to a decrease in Qm as ex- 
pected, because this is the direction of increase in ISW, 
which is needed given that we find the measured ISW 
exceeds WMAP3 prediction. The effect is smaller than 
expected because moving along the WMAP degeneracy 
line in the direction of decrease in ilm also requires an 
increase in h and decrease in erg, both of which reduce 
ISW (see also Fig. [151). 

For ACDM + fix model, we improve significantly over 
what using CMB alone can do. In Fig. [16] we com- 
pare 1-D distributions of the and ilm when we use 
WMAP+ISW+WL versus using WMAP alone. The 
ISW effect, as discussed above, can constrain the change 
of gravitational potential of the Universe as it depends 
linearly on the change of growth factor of the potential 
{D{a)/a). For example, in the closed Universe model we 
plotted in Fig[T5l D{a)/a increases as redshift decreases, 
while in the other two models, D{a)/a would decrease 
as redshift decreases. As dcfi/drj has a different sign for 
the closed Universe model on WMAP degeneracy curve 
as compared to the open and the fiat universe model on 
the same curve, the sign of the ISW effect changes too. 
In Fig. [15] we plot the predicted ISW signal using the 
low-redshift LRG and the high-redshift quasar distribu- 
tion for 3 different Universes along the WMAP degener- 
acy curve. As expected closed model differs drastically 
from open and fiat models. We also see that for LRG 
there is not much difference between fiat and open mod- 
els even though the latter has Q.ra = 0.088 compared to 
rim = 0.24, but the increase in ISW induced by D{a)/a 
is compensated by the reduction caused by other param- 
eters such as h and a^. The differences between the two 
are more significant for the high-z quasar redshift distri- 
bution. ISW effect breaks the WMAP degneracies be- 
tween ^Ik and f^A (or fim)- In Fig. [17] we show the 2-D 
contour plots of this set of parameters to demonstrate 
how our analysis improves the constraints on these pa- 
rameters. 

Finally, we look at the ACDM + w model where we 
look for better constraints on dark energy equation of 
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FIG. 16: ACDM + Oa' model: the 1-D distributions of 
and f2m. The solid (dashed) line represents constraints from 
using WMAP+ISW+WL (WMAP alone). 



TABLE IX: Comparing the constraints for several parameters 
with and without Weak Lensing of CMB in ACDM + Q.k 
cosmological model. The limits shown are mean and standard 
deviation for each of the parameter. 

Parameter Limits (CMB+ISW+WL) Limits (CMB+ISW) 
Qk -0.068 ±0.019 -0.0073 ± 0.020 

0.746 ± 0.059 0.745 ± 0.065 

n„ 0.261 ± 0.075 0.263 ± 0.083 
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FIG. 18: ACDM + w model: the 2-D distribution of 
Q.A and w (68% and 95% confidence contours shown). 
The solid (dashed) line represents constraints from using 
WMAP+ISW+WL (WMAP alone). 




K 

FIG. 17: ACDM + model: the 1-D distribution of 
and the 2-D distribution of and fix (68% and 95% confi- 
dence contours shown). The solid (dot-dashed) line represents 
constraints from using WMAP+ISW+WL (WMAP alone). 



state (w). The constraint on w is modestly improved, 
since the dark energy equation of state changes the 
growth factor along the WMAP degeneracy curve, thus 
the evolution of the gravitational potential. We also see 
that there is a tilt of flA towards lower value when we 
combine WMAP with ISW and WL effects. We also plot 
the 2-D contours for the flA and w in FigfTSl 

As mentioned above, WMAP 3-year model predicts 
ISW amplitude that is about two sigma below our con- 
straints and this is also the case for the best fit ACDM 
model (which is almost the same as WMAP 3- year). 
Adding curvature or dark energy equation of state does 
not reduce this discrepancy either and in both cases these 
two parameters are not needed to improve the fit. While 
reducing matter density goes in the desired direction of 
increasing ISW in cross-correlations, such models also in- 
crease the CMB power at large scales through the ISW 
auto-correlation power, which is in disagreement with the 
low power observed on large scales in WMAP. For exam- 
ple, we find that there are models with Qm = 0.18 which 
improve the fit to ISW data by 13 relative to the best 
fit ACDM + fix model, but at the same time make the 
WMAP fit worse by 15. There is thus some mild 
tension between low power in WMAP at low I and the 
high ISW power we measure, but it is a tension that can- 
not be removed by simple extensions of parameter space 
explored here. As this is only a two sigma effect there 
is a considerable probability that it is just a statistical 
fiuctuation. 
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TABLE X: The percentiles of the posterior distribution (2.5%, 16%, 50%, 84% and 97.5%) on cosmological parameter for each 
model with the CMB only ("C") and also including the ISW and weak lensing likelihood functions ("I"). Ho is in kms~^ Mpc~^. 
For a Gaussian distribution these percentiles correspond approximately to —2a, —Icr, central, +la, and +2a values. Note that 
for the 7-parameter chains with CMB only there are significant prior effects in the CMB degeneracy direction. 



Parameter C(2.5%) 


C(16%) C(50%) 


C(84%) 


C(97.5%) 


1(2.5%) 


1(16%) 


1(50%) 


1(84%) 1(97.5%) 


ACDM, 6 parameters 




0.0208 


0.0214 


0.0222 


0.0229 


0.0236 


0.0208 


0.0215 


0.0222 


0.0229 


0.0236 




0.0901 


0.0976 


0.105 


0.113 


0.121 


0.0901 


0.0970 


0.104 


0.111 


0.119 


T 


0.0312 


0.0612 


0.0911 


0.121 


0.151 


0.0359 


0.0662 


0.0956 


0.125 


0.154 


ris 


0.929 


0.943 


0.959 


0.976 


0.993 


0.929 


0.944 


0.960 


0.977 


0.994 




0.684 


0.724 


0.760 


0.793 


0.822 


0.698 


0.734 


0.766 


0.796 


0.822 




0.178 


0.207 


0.240 


0.276 


0.316 


0.178 


0.204 


0.234 


0.266 


0.302 




0.670 


0.717 


0.767 


0.816 


0.863 


0.671 


0.715 


0.763 


0.810 


0.855 


Ho 


67.0 


69.9 


72.9 


76.3 


79.7 


67.9 


70.6 


73.5 


76.6 


79.8 








ACDM + 


Oa', 7 parameters 












0.0203 


0.0211 


0.0218 


0.0226 


0.0233 


0.0206 


0.0213 


0.0221 


0.0229 


0.0236 




0.0916 


0.0990 


0.107 


0.115 


0.123 


0.0900 


0.0968 


0.104 


0.112 


0.120 


T 


0.0269 


0.0546 


0.0836 


0.113 


0.142 


0.0330 


0.0637 


0.0934 


0.123 


0.152 




-0.147 


-0.115 


-0.0499 


-0.00574 


0.0150 


-0.0515 


-0.0235 


-0.00395 


0.0103 


0.0201 




0.917 


0.932 


0.948 


0.966 


0.984 


0.925 


0.941 


0.958 


0.976 


0.993 




0.332 


0.437 


0.606 


0.745 


0.821 


0.610 


0.691 


0.754 


0.802 


0.837 




0.166 


0.262 


0.445 


0.678 


0.804 


0.148 


0.190 


0.250 


0.330 


0.436 


0-8 


0.648 


0.690 


0.738 


0.788 


0.839 


0.663 


0.709 


0.758 


0.807 


0.857 


Ho 


40.5 


43.6 


53.8 


69.5 


86.6 


54.0 


62.1 


71.0 


81.3 


92.0 










ACDM -+ 


w, 7 parameters 












0.0208 


0.0215 


0.0222 


0.0231 


0.0239 


0.0207 


0.0214 


0.0222 


0.0230 


0.237 




0.0900 


0.0981 


0.106 


0.114 


0.122 


0.0906 


0.0975 


0.105 


0.112 


0.120 


T 


0.0294 


0.0600 


0.0894 


0.119 


0.149 


0.0347 


0.0647 


0.0940 


0.123 


0.153 


W 


-1.731 


-1.457 


-1.031 


-0.573 


-0.240 


-1.646 


-1.401 


-1.006 


-0.704 


-0.425 


ris 


0.927 


0.943 


0.960 


0.981 


1.010 


0.928 


0.943 


0.960 


0.978 


0.998 


Oa 


0.457 


0.617 


0.764 


0.844 


0.870 


0.546 


0.672 


0.778 


0.845 


0.871 




0.130 


0.156 


0.235 


0.383 


0.543 


0.128 


0.155 


0.220 


0.328 


0.454 


crs 


0.437 


0.613 


0.776 


0.919 


1.032 


0.540 


0.659 


0.781 


0.898 


1.00 


Ho 


47.8 


57.9 


73.8 


90.7 


98.5 


53.2 


62.4 


75.6 


90.2 


98.2 



VIII. DISCUSSION 

The main goal of this paper is to perform a full anal- 
ysis of the integrated Sachs- Wolfe effect using the cross- 
correlations between WMAP CMB maps and maps of 
large scale structure. In contrast to previous work on 
this subject we place less emphasis on establishing a de- 
tection of ISW and more emphasis on developing a tool 
with which cosmological models can be compared to the 
data in a close to optimal fashion. For this reason we 
only select the data sets that can be reliably used to- 
wards this goal, as discussed in more detail below. The 
redshift range of the datasets we use is between and 
2.5. We use optimal weighting of the data both in angu- 
lar space and in redshift space to extract the maximum 
amount of information, taking into account properly the 
correlations between them. Our final product is the like- 
lihood function with which different cosmological models 
can be compared to each other. 

As the ISW effect is both a probe of cosmologi- 
cal parameters and a consistency test of the standard 
ACDM cosmology, there have been significant previous 
efforts made to detect it. A number of different analy- 



sis methods have been used and the WMAP data have 
been cross-correlated with several samples. These in- 
clude the 2MASS XSC; several SDSS samples including 
magnitude-Hmited galaxy samples, LRGs, and quasars; 
the NVSS; and the HEAO hard X-ray map. Most of 
these samples (or samples with similar spatial coverage 
and redshift range) are included in the present work, but 
not all. Here we compare our analysis with the previous 
work and comment on the reasons for our choice of data 
sets. 

1. Near-infrared galaxies (2MASS). The 2MASS 
galaxies are useful for ISW due to high sky cov- 
erage and the ability to see closer to the Galactic 
plane in the near-IR than in the optical. However 
they can only probe the lowest redshifts {z < 0.2). 
Afshordi et al. and Rassat et al. have mea- 
sured the ISW signal using the 2MASS sample and 
we dehberately cut our 2MASS sample into bright- 
ness bins such as theirs so that we can compare 
the results. We find that our measured signal from 
2MASS is very similar. We do however derive cos- 
mological constraints using a Markov chain (which 
fits all the cosmological parameters instead of just 
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JIa) from these samples. We also take into account 
(albeit in a crude way) the redshift dependence of 
the bias resulting from seeing all nearby galaxies 
but only the brightest and most biased galaxies at 
z > 0.1. 

Optical galaxies (SDSS, APM). Wide-angle multi- 
color galaxy surveys such as SDSS open almost lim- 
itless possibilities for constructing galaxy samples, 
and many of these samples have been used in pre- 
vious ISW work. Most work so far has been on 
either flux-Hmited samples [H, US], which have a 
broad redshift distribution, or photometric LRGs 
[l^ . [isl . [isl . [20| . which can be seen to larger dis- 
tances and for which it is easier to construct re- 
liable photo-z cuts. In SDSS, photometric LRG 
samples oversample the Hnear density field in the 
redshift range 0.2 < z < 0.6 and the lower redshifts 
are covered by 2MASS, so the flux- limited galaxy 
samples would be redundant in terms of volume 
for our study; we therefore did not include them. 
Our LRG samples cover the largest solid angle to 
date of any SDSS ISW analysis (6641 deg^) and 
for the purposes of cosmological analysis are split 
into two photo-z slices. Fosalba & Gaztafiaga 
have also used galaxies from the Automated Plate 
Measuring (APM) survey which adds ~ 4300 
deg^ in the Southern Hemisphere inaccessible to 
SDSS. Their APM sample has a typical redshift 
z « 0.15 and thus would add some information be- 
yond the most distant of our 2MASS samples. Con- 
sidering that APM area is only 16% of 2MASS and 
that it only marginally extends the redshift range 
we have not used APM in our analysis. However 
adding a deeper galaxy survey in the South, com- 
parable to or deeper than SDSS, would be valu- 
able for improving ISW constraints. Overall sig- 
nal to noise from SDSS LRG galaxies is about 3<j, 
most of which comes from the higher redshift sam- 
ple around z ~ 0.5. 

Optical quasars (SDSS). Photometrically selected 
quasars can probe large-scale structure at much 
higher redshifts than "normal" galaxies because 
they are bright enough to be seen in wide-angle 
surveys (such as SDSS) even at z ^ 2. The only 
ISW analysis with quasars so far has been that 
of Giannantonio et al. [2l|, who cross-correlated 
WMAP with a sample of photometric quasars from 
the SDSS. Our analysis uses similar selection cri- 
teria, but we have used photo-z cuts to eliminate 
most of the lower-redshift objects, and used a com- 
bination of spectroscopic data and angular cluster- 
ing to constrain b * dN/dz taking into account the 
multimodal nature of the photo-z failures. We also 
slice our quasars into two photo-z bins. Despite 
these improvements we find that the significance is 
only 1.3(7 (1.24(7 for the high redshift sample with 
z > 1), and we therefore do not confirm that the 



2.1(7 signal seen in [2l| comes from z > 1. 

4. Radio sources (NVSS). There have been several 
past WMAPxNVSS ISW analyses [H, [13, Sill, 
taking advantage of the high redshift (compared 
to most optical samples) and wide sky coverage of 
the NVSS. We have used the angular power spec- 
trum whereas the previous works have used corre- 
lation functions or wavelet coefficients. However, 
the most important difference between our analy- 
sis and the previous result is that we fit 6 * dN/ dz 
from cross-correlations rather than using the Dun- 
lop & Peacock model (ttI I for the redshift distribu- 
tion and assuming constant bias. This is important 
as we find the fit 6 * dN/ dz looks very different (see 
Fig. [T3|) . All of these studies, including ours, have 
found positive cross correlations at the ~ 3(7 level. 
However, the interpretation of this result depends 
sensitively on one's ability to measure b * dN/dz 
and this is where we believe our analysis is an im- 
provement upon previous efforts. 

5. Hard X -ray background (HE AO). Boughn & Crit- 
tenden [15| have used the HEAO hard X-ray map 
(90| for ISW cross-correlation. The background is 
due mainly to unresolved (by HEAO) active galac- 
tic nuclei and hence traces large scale structure 
at redshifts of order unity. This, combined with 
the all-sky nature of HEAO, is beneficial for ISW 
projects. However, we decided not to add in HEAO 
sample to our analysis for several reasons. First 
is the difHculty in understanding the 6(z) * dN/dz 
of the sample (we use the general notation dN / dz 
here even though for unresolved X-ray fiux it would 
be more accurate to write dF/dz). Only ^ 75% of 
the background is resolved b y C handra into sources 
with measured redshifts [9ll.l93|. and we have little 
guidance on where to place the other 25%. Even 
if we knew dN/dz perfectly, this does not tell us 
b * dN/dz: the modeled dN/dz spans the range 
< z < 3 and it is unlikely that the bias would 
be even approximately constant over this range. 
An alternative is to fit for their bias and redshift 
distributions up to high z using a cross-correlation 
method similar to that done for NVSS in Sec. IVl 
Unfortunately HEAO has FWHM of - 2° and does 
not resolve individual sources, so we would have to 
fit the data to the model without small-scale in- 
formation, which loses signal-to-noise on the cross- 
correlation very rapidly. A secondary reason is that 
there is considerable overlap between HEAO and 
NVSS, so it is likely that the two trace partly the 
same structure, and thus the improvement in ISW 
constraint is not as large as adding two indepen- 
dent data sets. We note that it may make sense 
to include the hard X-ray maps in parameter esti- 
mation in the future if a robust determination of 
b * dN/dz becomes available. 

In summary, we believe we used most of the available 
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large scale structure data useful for ISW analysis. This 
not only updates previous work on ISW effect [ill, [13, 

[B, Q [il, M, [13, in, E H [U, [13, [H, [H, [11 , but is 

also the first one that attempts to do the tomography of 
ISW, in the sense of encompassing a wide redshift range 
via our mass tracers going from the local Universe to 
z ^ 2.5, while reducing the amount of overlap in area and 
redshift as much as possible. We have argued that many 
of the previous measurements have a considerable overlap 
in redshift and area, which means that they cannot be 
combined independently and that the effective redshift of 
the sample is not necessarily the redshift from where most 
of the ISW signal is coming from. Our analysis, while 
attempting to minimize the overlap in the first place, 
takes the residual correlations into account explicitly via 
the construction of the full covariance matrix. We note 
that Giannantonio et al. [9^ are also pursuing an ISW 
tomography analysis, with somewhat different choices of 
LSS samples and cross-correlation methodologies. 

We spend a significant fraction of our analysis obtain- 
ing the correct redshift distributions for all of the sam- 
ples. To be more accurate, it is the b * dN/dz that we 
constrain for all samples. The signature of ISW effect is 
highly affected by the redshift distribution of the tracer, 
and thus one would need to have an accurate idea of what 
the redshift distribution is in order to interpret the corre- 
lation. Apart from employing spectroscopic datasets that 
overlap in magnitude range and sky coverage, we corre- 
late the tracer samples with one another so as to obtain 
the b * dN/dz for some of the samples. This is mainly 
possible because LRGs have relatively good photometric 
redshifts and so we correlate the LRGs with other over- 
lapping datasets to determine what are the b * dN/dz 
at the redshift range that LRGs cover. In addition, we 
account for redshift-dependent bias in 2MASS and for 
the multimodal error distributions for the quasars. We 
also made the first determination of 6 * dN/dz for NVSS 
sample which is not based simply on a theoretical model 
fitting the luminosity function. 

Correlations of mass tracers with the CMB sky can be 
caused not only by the ISW effect, but also by other cos- 
mological effect such as thermal SZ, Galactic foregrounds 
and extinction, and extragalactic point sources. We pro- 
vide an estimate for all these effects and only include 
the scales deemed reliable, where the contamination is 
subdominant or negligible. 

We report a detection of 3.7(7 of the ISW effect com- 
bining 2MASS, SDSS, and NVSS with WMAP data. We 
make a joint analysis of all samples by constructing a reli- 
able covariance matrix including cross-correlations of dif- 
ferent samples, which is needed for cosmological param- 
eter fitting. We combine our ISW correlation functions 
with weak lensing of the CMB (Paper II) to derive cosmo- 
logical constraints on three different cosmological models: 
(i) the "vanilla" ACDM model, (ii) ACDM+Hk, and (iii) 
ACDM+ui. We find a slight improvement of our mea- 
surement of w in model (iii) over the measurement made 
by CMB alone: w = -l.Olt^j^ instead of -1.03l^;4^. 



The constraining power of our analysis is however most 
prominent in determining that curvature of the Universe: 
for CMB+ISW+WL we find = -0.004^°:°^;^ instead 
of -O.OSOto Jjgs for CMB alone. These constraints are 
not as tight as that obtained by some other methods, 
such as combining the CMB with baryonic oscillations 
or with supernovae (gsI. [9^. [93| . but it is subject to very 
different systematics. It is thus reassuring that all of 
them are consistent with each other. Even more impor- 
tantly, there are other models where ISW can be cru- 
cial in distinguishing them from standard ACDM, such 
as f{R) models in which the growth of structure is not 
fixed by the background geometry [l^l- Some of these 
models may already be inconsistent with our ISW signal; 
we plan to present such constraints in a future paper. 
These constraints should improve further in the future 
with deeper galaxy surveys that should reach the cosmic 
variance limit out to z ~ 1 — 2, and future CMB data 
that enables lower-noise lensing reconstruction. 

Finally, we would like to note that we plan to release 
a pacl-cage for calculating ISW likelihood function given 
the datasets and cosmological parameters. This will be 
described further in the documentation for the package. 
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APPENDIX A: NVSS REDSHIFT WINDOW 
FUNCTIONS 



In Sec. IV D 3| we fit a constant /nvss to the cross- 
spectra of NVSS and the other LSS samples. We argued 
that this procedure gave an estimator /nvss whose ex- 
pectation value was given by Eq. l(34)) . The purpose of 
this appendix is to prove this equation and construct the 
functional form yV(z). 

We begin by noting that we have measured cross- 
spectra Cl'^^^^ and their covariance matrix Y.^,. The 
theoretical cross-spectrum is on the other hand simply 
the Limber result. 



xP k 



e + 1/2 



dz 



r{z)YH{z) 
hYSsiz)K.liz)dz, (Al) 



where r{z) is the comoving angular diameter distance 
and kI{z) is defined by the equivalence in the second 
line. The fitting procedure for the constant /nvss is 
to minimize 

= ^p-l]«<((?r^^^-/NVSsi^^) 



^/Ai,NVSS r r^i ^ 

X(^li -jNVSS^«'j, 



(A2) 



cross-spectra. The minimum value of is obtained for 



/i 



NVSS 



(A3) 



Since the (J^^^^^^^ have expectation value given by 



Eq. IjAip . we have 



(/i 



NVSS/ 



(A4) 

This proves Eq. (jM)) and shows that the window function 
is 



W{z) = 



j:,,,[E^-^WKl,^,l{z[ 



(A5) 



APPENDIX B: ERROR BARS ON GALAXY-CMB 
CORRELATIONS 



where Ki = K\{z)dz and C\' are the measured 



The purpose of this appendix is to discuss our choice 
of the Monte Carlo ("MCI") estimator for the error bars 
on the estimator, and then give a crude estimate 
for the possible biases that are induced by its use. As 
mentioned in the main text there are two biases: the 
correlation bias (because the galaxies and CMB are cor- 
related and MCI does not take this into account) and 
a realization bias (since we have only one realization of 
the galaxies) . The correlation bias is deterministic in the 
sense that the error bar is always underestimated in ev- 
ery ^-bin. The realization bias is more subtle: liCf^ = 0, 
then the MCI estimator returns an unbiased estimate of 
a'^{Cf^). However the £-bins where the error is under- 
estimated are weighted more heavily than those where it 
is overestimated, resulting in a final error bar on cosmo- 
logical parameters that is biased low. 

We consider each of these issues separately in a toy 
model. The toy model has the following assumptions: 

1. We are computing cross-spectra Cf^ in M ^-bins 
(call these cross-spectra xi...xn)- 

2. The galaxies and CMB temperature are Gaussian 
random fields. (We are at Hnear scales where large 
scale structure is Gaussian; the systematics may 
not be.) 

3. In the ith ^-bin, there are Ni galaxy modes, and 
all of the CMB modes in this region are observed. 
(This is true except that NVSS goes slightly closer 
to the Galactic plane than WMAP.) We ignore 
mode coupling at the boundaries, i.e. each mode is 
ascribed to a single ^-bin. 

4. We are fitting the cross-correlation data to some 
parameter, say an amplitude A of some template 
ti. More generally, when the ISW effect is essen- 
tially constraining one direction in parameter space 
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with all others constrained by the CMB alone (the 
case with ACBM+flK and ACDM+w models here) 
the template would be dxi/dp where p parameter- 
izes the CMB-degenerate direction. The fit is done 
using the Monte Carlo covariance matrix. 

5. The objective is to determine what is the ratio of 
the "true" error bar on A to that derived from the 
fitting procedure. 

Within these assumptions, we evaluate the correlation 
bias i?i and realization bias i?2, which we define to be 
the ratio of true to estimated variance. We find, us- 
ing correlation coefficients and numbers of modes for the 
worst-case bins, that Ri ~ 1.02 and {R2) = 1.11. This 
corresponds to ~ 6% underestimation of the error bars 
in the worst case, which is negligible. 



where the correlation coefficient is 

p, = , ^' (B6) 

For the fiducial cosmology and the bins that we used, 
the maximum predicted correlation coefficient is 0.067 
(LRGl, £ = 18). This would suggest an underestimate of 
the error bar by a factor of i?i = 1.0044. For some cosmo- 
logical models, such as those with lower Qm or higher cts, 
the correlation coefficient could be larger. Indeed there is 
some evidence for this: we observe an overall ISW ampli- 
tude of 2.2 ± 0.6 times the prediction. If we multiply the 
correlation coefficient p by 2.2 then the underestimate of 
the error bar grows to Ri = 1.02; even this is negligible. 



1. Correlation bias 

We will introduce the notation C^^ — C^^ +n~^ for the 
galaxy power spectrum including Poisson noise, and for a 
matrix Cov we will write Cov~-^ to mean the ij element of 
Cov~^ rather than the reciprocal of Covy . (In the cases 
considered in this appendix the covariance matrices are 
diagonal so this distinction will not matter.) We will 
also use the shorthand Cf for the galaxy-temperature 
cross-spectrum in the ith bin. 

The estimator for the cross correlation is 



(Bl) 



where a = l...Ni is a mode index. 

The true uncertainty in Gaussian theory, using inde- 
pendence of modes, is 



2. Realization bias 

Having taken into account the correlation bias, we now 
consider the case where the cross-correlation coefficient 
is small (p ^ 1). In this case, the covariance matrix of 
the that we obtain from the CMB Monte Carlos is 
unbiased. The realization bias comes from the fact that 
we invert the covariance matrix, and unbiased Cov does 
not imply unbiased Cov~l. 

The true covariance matrix of the estimator Eq. (|Bip 
for Xi is 



Covi 



Cov(a;,,a::_,) 



(B7) 



where Cf^ and Cj"^ are the true (ensemble-averaged) 
galaxy and CMB power spectra, including Poisson noise 
for the galaxies. The estimated covariance matrix is in- 
stead given by Eq. (|B4p . Now define the number 



Cov," 



Cov(a;i, Xj] 



"f[c!^cr+{cfn (B2) 



However the MCI procedure gives 

CoVij = ^ gagg(rar/3)MontcCarlo, (B3) 

where the a modes are in bin i, the beta modes are in 
bin j, and g is the actual realization of the galaxies. Sim- 
plifying with CMB covariance matrix gives 



(B4) 



Note that this is diagonal, even though we have only 
Monte-Carloed one of the data sets. 

In the presence of a nonzero cross-correlation, the MCI 
covariance matrix is biased: 



Gov,;, 



= 1 



p1, 



(B5) 



Gov,;. 

Govi. 



(B8) 



which is the ratio of the estimated to true variance in a 
given bin. This simplifies to 



(B9) 



i.e. it is a distribution with Ni degrees of freedom, 
divided by the number of degrees of freedom. In partic- 
ular {uij = 1: the covariance matrix is unbiased, but we 
have from distribution theory {y~^) = (1 — 2/N)~^ 
and {y-'^) = (1 - 2/7V)-i(l - 4/iV)-i. 

However what we really want to compare are the true 
and estimated errors on the parameter A. The estimate 
A of the amplitude A is 



^ ^ Gov^j tjx^ 
Govy titj 



(BIO) 
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Its estimated variance is 



1 



Var(A) = _^ 

Its true variance is 



CoVjj titj 



(Bll) 



Var(i) 



thCov^^ Covy Cov J. tk 



Covy Utj 



(B12) 



(Note that Xi and Covy are uncorrelated because the 



probabiUty distribution is symmetric in Tq 



-Tq, un- 



der which Xi changes sign but Covy does not.) The ratio 
is 

^ _ Var(i) thCovJcovijCovj^tk 
-"-2 = ^r—, ^- — . ,_i • (i^ijj 



Var(^) 



Covjj titj 



Using the definition of yi and diagonahty of the matrices 
Gov and Gov, 



l„.-2 



i?2 



(B14) 



We now consider two Hmiting cases. If we have a single 
i'-bin, then the ratio is i? = yf^ and 



{R2) = 



1 - 2/N, 



(B15) 



If we have many £-bins contributing then the sums go to 
their mean values and we get 



{R2 



1 



1- w 



(B16) 



if there were the same number of modes in each £-bin. 
This is larger than Eq. I|B15|) because with only a single 
bin there is then no possibility for the amplitude estima- 
tor to re-weight the bins to take advantage of the ones 
with smaller estimated variance. 

The number of modes per bin is, in the limit of negli- 
gible mode coupling, 



Ni — [(^max + 1) ^min]/f 



'sky 



(B17) 



This is 40 for the lowest-^ 2MASS bin that we use in pa- 
rameter fits, 77 for the lowest-^ LRGs, 70 for the quasars, 
and 40 for NVSS. To be pessimistic, if we take Eq. (|BT61) 
for all these cases the worst number we get is {Ra) — 1-11, 
which means that in this pessimistic case we have under- 
estimated the error bar (a) on the cross-correlation by 
5%. In reality much of the constraint comes from higher- 
£ bins where Ni is greater, so this should be taken as an 
upper limit. 



APPENDIX C: SIGNAL TO NOISE ESTIMATE 
AND UPPER LIMIT ON FOREGROUND 
CONTAMINATION 



To assess the statistical signal to noise we look at cor- 
relation between the galaxy overdensity and the tem- 
perature of the GMB (C'^of,^), which is the data vec- 
tor called d. We also need its inverse covariance ma- 
trix (C~^) and the theoretical prediction, which we can 
model as amplitude A times a fiducial model t. To assess 
the possible contamination from foregrounds, tSZ, point 
sources etc. to our signal we need to estimate the as- 
sociated cross-correlation contamination (/). For exam- 
ple, for foregrounds we look at the correlation between 
the galaxy overdensity and the foreground temperature 
cf '^^^ (which is what we calculated using models such 
as described by Eq. [44|) . 

Gonsider the usual analysis, where we are trying to 
fit for A given d, t and C~^: 



= {d- At) ■ C-\d- At). 
We minimize cind get 



A = 



d-C-H 
t- c-if 



(CI) 



(G2) 



This is Eq. l(45|) and the associated variance is given by 
Eq. l(46|) . The ratio of estimated amplitude to its variance 
is the estimated signal to noise. 

Since the total signal is a sum of the true signal and 
contamination such as foreground, tSZ or point sources, 
the latter contribute to the signal to noise. 



f-c-H f-c-^t 



AA 



Vt- C-if Vf • c-if 



/■c-i/ 



(G3) 



While we could use this expression to estimate the pos- 
sible contamination we can make it less dependent on 
the weighting by theoretical model by using the Gauchy 
inequality, here written in the (primed) diagonal ba- 
sis with eigenvectors normaHzed to eigenvalue, f ■ t' < 

\/t' -t'Jf ■ /', to derive from Eq. ([C3l) : 



(G4) 



We use this expression in our estimates of contamination; 
it represents an upper limit on the number of sigmas of 
contamination introduced by the foreground /. 
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